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Chapter 



Introduction 



From the seminal ideas of Feynman [1] and until now, quantum information and computation [2] 
has been a rapidly evolving field. While at the beginning, physicists looked at quantum mechanics 
as a theoretical framework to describe the fundamental processes that take place in Nature, it was 
during the 80's and 90's that people began to think about the intrinsic quantum behavior of our 
world as a tool to eventually develop powerful information technologies. As Landauer pointed 
out [3], information is physical, so it should not look strange to try to bring together quantum 
mechanics and information theory. Indeed, it was soon realized that it is possible to use the laws 
of quantum physics to perform tasks which are unconceivable within the framework of classical 
physics. For instance, the discovery of quantum teleportation [4], superdense coding [5], quantum 
cryptography [6,7], Shor's factorization algorithm [8] or Graver's searching algorithm [9], are some 
of the remarkable achievements that have attracted the attention of many people, both scientists and 
non-scientists. This settles down quantum information as a genuine interdisciplinary field, bringing 
together researchers from different branches of physics, mathematics and engineering. 

While until recently it was mostly quantum information science that benefited from other fields, 
today the tools developed within its framework can be used to study problems of different areas, like 
quantum many-body physics or quantum field theory. The basic reason behind that is the fact that 
quantum information develops a detailed study of quantum correlations, or quantum entanglement. 
Any physical system described by the laws of quantum mechanics can then be considered from the 
perspective of quantum information by means of entanglement theory. 

It is the purpose of this introduction to give some elementary background about basic concepts of 
quantum information and computation, together with its possible relation to other fields of physics, 
like quantum many-body physics. We begin by considering the definition of a qubit, and move then 
towards the definition of entanglement and the convertibility properties of pure states by introducing 
majorization and the von Neumann entropy. Then, we consider the notions of quantum circuit and 
quantum adiabatic algorithm, and move towards what is typically understood by a quantum phase 
transition, briefly sketching how this relates to renormalization and conformal field theory. We also 
comment briefly on some possible experimental implementations of quantum computers. 
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Chapter 0. Introduction 



What is a "qubit"? 

A qubit is a quantum two-level system, that is, a physical system described in terms of a Hilbert 
space C 2 . You can think of it as a spin-^ particle, an atom in which we only consider two energy 
levels, a photon with two possible orthogonal polarizations, or a "dead or alive" Schrodinger's 
cat. Mathematically, a possible orthonormal basis for this Hilbert space is denoted by the two 
orthonormal vectors |0> and |1). This notation is analogous to the one used for a classical bit, which 
can be in the two "states" or 1. Notice, however, that the laws of quantum mechanics allow a qubit 
to physically exist in any linear combination of the states |0> and |1). That is, the generic state \if/) 
of a qubit is given by 

|^> = o|0>+^|l>, (1) 

where a and B are complex numbers such that \a\ 2 + \B\ 2 = 1. Given this normalization condition, 
the above state can always be written as 

|«A> = ^> (cos (I) 10) + ^ sin (|) H>) . ( 2 ) 

where y, 9 and <f> are some real parameters. Since the global phase e iy has no observable effects, the 
physical state of a qubit is always parameterized in terms of two real numbers 8 and 0, that is, 

|^> = cos (I) |0> + sin (|) Il> ■ (3) 

The angles 9 and <p define a point on a sphere that is usually referred to as the Block sphere. Gen- 
erally speaking, it is possible to extend the definition of qubits and define the so-called qudits, by 
means of quantum <i-level systems. 

What is "entanglement"? 

The definition of entanglement varies depending on whether we consider only pure states or the 
general set of mixed states. Only for pure states, we say that a given state \if/) of n parties is entangled 
if it is not a tensor product of individual states for each one of the parties, that is, 

W> * |V1>1 ® |V 2 >2® •••® \v n ) n . (4) 

For instance, in the case of 2 qubits A and B (sometimes called "Alice" and "Bob") the quantum 
state 

l<A + > = ^(|0>a®|0> b + |1> a ®|1> b ) (5) 
V2 

is entangled since + \va)a ® \vb)b- On the contrary, the state 

\<f>) = \ (|0>a ® |0> B + |1>a ® |0> B + |0>a ®\\)b + IDa ® |1)b) (6) 



is not entangled, since 



|0> = |-^(|O>a + |1>a))®(-^(|O> b + |1> b )| . (7) 
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A pure state like the one from Eq|5]is called a maximally entangled state of two qubits, or a Bell 
pair, whereas a pure state like the one from EqQis called separable. 

In the general case of mixed states, we say that a given state p of n parties is entangled if it is 
not a probabilistic sum of tensor products of individual states for each one of the parties, that is, 



with {pk} being some probability distribution. Otherwise, the mixed state is called separable. 

The essence of the above definition of entanglement relies on the fact that entangled states of 
n parties cannot be prepared by acting locally on each one of the parties, together with classical 
communication (telephone calls, e-mails, postcards...) among them. This set of operations is often 
referred to as "local operations and classical communication", or LOCC. If the actions performed 
on each party are probabilistic, as is for instance the case in which one of the parties draws a ran- 
dom variable according to some probability distribution, the set of operations is called "stochastic 
local operations and classical communication", or SLOCC. Entanglement is, therefore, a genuine 
quantum-mechanical feature which does not exist in the classical world. It carries non-local correl- 
ations between the different parties in such a way that they cannot be described classically, hence, 
these correlations are quantum correlations. 

The study of the structure and properties of entangled states constitutes what is known as entan- 
glement theory. In this thesis, we shall always restrict ourselves to the entanglement that appears 
in pure states. We also wish to remark that the notation for the tensor product of pure states can 
be different depending on the textbook, in such a way that \va)a ® \vb)b = \va)a\vb)b = \va,vb)- 
An introduction to entanglement theory, both for pure and mixed states, can be found for instance 
in [10]. 

Majorization and the von Neumann entropy 

Majorization theory is a part of statistics that studies the notion of order in probability distributions 
[11-14]. Namely, majorization states that given two probability vectors x and y, the probability 
distribution y majorizes x, written as x < y, if and only if 



where {p^} is a set of probabilities and {P^) is a set of permutation matrices. The above definition 
implies that the probability distribution x is more disordered than the probability distribution y, since 
it can be obtained by a probabilistic sum of permutations of y. More details on majorization theory, 
which is often used in this thesis, are given in Appendix A. 

Majorization theory has important applications in quantum information science. One of them 
is that it provides a criteria for the interconvertibility of bipartite pure states under LOCC. More 
concretely, given two bipartite states \^ab) and \4>ab) for parties A and B, and given the spectrums 

and p0 of their respective reduced density matrices describing any of the two parties, the state 
\^jab) may be transformed to \4>ab) by LOCC if and only if [15] 




(8) 



k 




(9) 



k 



P>f> <P<I> ■ 



(10) 
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Chapter 0. Introduction 



An important theorem from classical information theory that plays a role in the study of en- 
tanglement is the so-called theorem of typical sequences. In order to introduce it, let us previously 
sketch some definitions. Consider a source of letters x which are produced with some probability 
p(x). The Shannon entropy associated to this source is defined as H - - J^ x P( x ) l°g2 PC*)- Given a 
set of n independent sources, we say that a string of symbols {x\, x 2 , ...,*„) is e-typical if 

2~n(H-e) < p(xi , X2 , . . . , X J < 2 -'<^> , (11) 

where p(x\, x 2 , ■ ■ ■ , x n ) = p{x\)p{x2) ■ ■ ■ p(x n ) is the probability of the string. The set of the e-typical 
sequences of length n is denoted as T(n, e). We are now in position of considering the theorem of 
typical sequences, which is composed of three parts: 

Theorem 0.1 (of typical sequences): 

• Given e > 0, for any 6 > and sufficiently large n, the probability that a sequence is e-typical 
is at least 1-5. 

• For any fixed e > and 5 > 0, and sufficiently large n, the number \T(n, e)| of e-typical 
sequences satisfies 

(1 - 5)2 n{H - e) < \T(n, e)| < 2 n(H+£) . (12) 

• Let S (n) be a collection of size at most 2 nR , of length n sequences from the source, where 
R < H is fixed. Then, for any 5 > and for sufficiently large n, 

^ p(xi,x 2 ,...,x„) < 8. (13) 

(x\,X2,—,x„) e S(n) 



It is not our purpose here to provide a detailed proof of this theorem (the interested reader is 
addressed for instance to [2]). We shall, however, make use of it in what follows. 

Let us introduce at this point a quantity which is to play a major role all along this thesis. 
Given a bipartite pure quantum state \4/ab), with reduced density matrices pa = ^b{Wab)(4'ab\) and 
Pb - trA{\4 f AB)('pAB\), the von Neumann entropy of this bipartition is defined as 

S =S (p A ) = -tr(p A log 2 p A ) = S (p B ) = -tr(p B log 2 Ps) , (14) 

where the equality follows from the fact that pa and ps share the same spectrum. This entropy is 
also called entanglement entropy, since it provides a measure of the bipartite entanglement present 
in pure states. To be precise, the entanglement entropy measures the optimal rate at which it is 
possible to distill Bell pairs by LOCC in the limit of having an infinite number of copies of the 
bipartite system. 

Let us explain how the above consideration works. Given the bipartite pure state \i^ab), we write 
it in terms of the so-called Schmidt decomposition: 



\<Pab) - ^ ^Ip(x)\xa)a\xb)b , 

X 



(15) 
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where the square p(x) of the Schmidt coefficients define the probability distribution that appears as 
the spectrum of the reduced density matrices for the two parties. The rc-fold tensor product \ip AB ) m 
can be written as 

I<Aab>®" = ^ ^p{x\)p{x 2 ) ■ ■ ■ p(x n )\x lA , x 2A , • • • , x nA ) A \x lB , x 2B , x nB ) B . (16) 

(xi ,x 2 ,-,x„) 

Let us now define a quantum state \<p n ) obtained by omitting in EqEDfhose strings (x\,x 2 , . . . , x n ) 
which are not e-typical: 

I0n> = ^ ^P(xi)p{x 2 ) ■ ■ ■ p(x n )\xi A , x 2A , x nA ) A \x XB , x 2B , x nB ) B . (17) 

(xi,X2,...,x n ) e T(n,e) 

Since the previous state is not properly-normalized, we define the state \cf>' n ) = \(p n )l ^(<f>n\4>n)- Be- 
cause of the first part of the theorem of typical sequences, the overlap between \ij/ AB ) m and \<p' n ) tends 
to 1 as n — > oo. Furthermore, by the second part of the theorem we have that \T(n, e)| < 2 n(H+ ^ = 
2 n< - s+£ \ Given these properties, a possible protocol to transform copies of the state \i// AB ) into Bell 
pairs by means of LOCC reads as follows: party A may convert the state \tp AB ) m into the state \(p' n ) 
with high probability by performing a local measurement into its e-typical subspace. The largest 
Schmidt coefficient of \<f> n ) is 2~ n(S ~ £ ^ 2 by definition of typical sequence, and since the theorem of 
typical sequences also tells us that 1 - 8 is a lower bound on the probability for a sequence to be 
e-typical, the largest Schmidt coefficient of \<p' n ) is at most 2 _ " (S_e ^ 2 / Vl - 5. Let us now choose an 
m such that 

2~n(S-e) 

< T m . (18) 

1 -5 

Then, the spectrum of the reduced density matrices for A and B are majorized by the probability 
vector (2" m , 2~ m , 2" m ) r , and therefore the state \<p' n ) can be transformed into m copies of a Bell 
state by means of local operations and classical communication. More specifically, in the limit 
n — > oo the ratio m/n between the number of distilled Bell pairs and the original number of states 
exactly coincides with the entanglement entropy S . 

It is possible to see that the above distillation protocol is optimal, that is, it is not possible to 
distill more than nS Bell pairs from a total of n copies of a bipartite pure state in the limit n — > oo. 
Because of this property, the von Neumann entropy is also called the distillable entanglement of a 
pure bipartite system. Furthermore, it is possible to see that the entropy S coincides as well with 
the entanglement of formation of bipartite pure states, which is the optimal ratio m/n describing the 
number m of Bell pairs that are required to create n copies of a given bipartite pure state by means 
of LOCC, in the limit n — » oo. The von Neumann entropy constitutes then a genuine measure of the 
bipartite entanglement that is present in a given pure quantum state. 

Quantum circuits and adiabatic quantum algorithms 

Much in analogy to the situation in classical computation, where it is possible to define a compu- 
tation by means of logic gates applied to bits, a quantum computation may be defined in terms of 
a set of unitary gates applied to qubits. These unitary gates may either be local, acting on a single 
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Figure 1: Quantum circuits representing the action of a Hadamard gate on a single qubit and a 
controlled-not gate on two qubits. The controlling qubit is denoted by a black dot, and the controlled 
qubit is denoted by the symbol ©. 

qubit, or non-local, acting on several qubits at a time. An important example of a local gate is given 
by the so-called Hadamard gate: 

U » = T^ -.)• < 19 > 



which acts on the two-dimensional Hilbert space of a single qubit such that 

U H \0) = -U|0> + |1» 
V2 

u H \i) = -L(|o>-|i». 

V2 

Also, an important example of a non-local gate is the controlled-not gate Ucnot- 



(20) 



Ucnot - 



(\ 0\ 

10 

1 

{0 1 



(21) 



acting on the four-dimensional Hilbert space of two qubits such that 



U C not\0,0) = |0,0) 

UcnotWA) = 10,1) 

UcnotIIO) = |1, 1> 

U CNO t\\,\) = |1,0). 



(22) 



In the example of the controlled-not gate, the first and second qubits are respectively called the 
controlling qubit and the controlled qubit, since the action of the gate on the second qubit depends 
on the value of the first one. It is possible to define more general controlled gates similarly to the 
controlled-not gate, namely, if the controlling qubit is in the state |0) nothing is done on the second 
one, whereas if the controlling qubit is in the state 1 1 ) then some local unitary gate acts on the second 
qubit. The application of the different unitary gates that define a quantum computation on a system 
of qubits can be represented in terms of quantum circuits, such as the ones from Fig^and Fig|2] In 
a quantum circuit each wire represents a qubit, and the time flows from left to right. 

Independently of quantum circuits, it is possible to define alternative models to perform quantum 
computations, such as the adiabatic model of quantum optimization [16]. The adiabatic quantum 
algorithm deals with the problem of finding the ground state of a physical system represented by 
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Figure 2: A possible quantum circuit of 5 qubits composed of Hadamard and controlled-not gates. 
Some measurements are performed on the qubits at the end of the quantum computation. 

its Hamiltonian Hp. The basic idea is to perform an interpolation in time between some easy-to- 
build Hamiltonian Hq and Hp, such that if the initial state of our system is a ground state of Hq, 
we may end up in a ground state of Hp with high probability after evolving for a certain amount 
of time, as long as some adiabaticity conditions are fulfilled. For example, we could consider the 
time-dependent Hamiltonian 

H(t) = (\-^Ho + jH P , (23) 

where t € [0, T] is the time parameter, T being some computational interpolation time. If g m! -„ 
represents the global minimum along the evolution of the energy gap between the ground state and 
the first excited state of the system, the adiabatic theorem implies that, if at t - the system is 
at ground state of Ho, in order to be at the ground state of Hp at time T with high probability it 
is required that T ~ 1 /g 2 . . The scaling properties with the size of the system of the minimum 
energy gap controls then the computational time of the quantum algorithm. Actually, the fact that 
the system evolves through a point of minimum gap implies that it approaches a quantum critical 
point, to be defined in what follows. A more detailed explanation of adiabatic quantum algorithms 
is given in Chapter 4. 

Quantum criticality in quantum many-body systems 

A quantum phase transition is a phase transition between different phases of matter at zero temper- 
ature. Contrary to classical (also called "thermal") phase transitions, quantum phase transitions are 
driven by the variation of some physical parameter, like a magnetic field. The transition describes 
an abrupt change in the properties of the ground state of the quantum system due to the effect of 
quantum fluctuations. The point in the space of parameters at which a quantum phase transition 
takes place is called the critical point, and separates quantum phases of different symmetry. 

Some properties of the system may display a characteristic behavior at a quantum critical point. 
For instance, the correlators in a quantum many-body system may decay to zero as a power-law 
at criticality, which implies a divergent correlation length and therefore scale-invariance, while de- 
caying exponentially at off-critical regimes. Since quantum correlations are typically maximum at 
the critical point, some entanglement measures may have a divergence. The ground-state energy 
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may display non-analyticities when approaching criticality, and the energy gap between the ground 
state and the first excited state of the system may close to zero. Our definition of quantum phase 
transition is very generic and does not necessarily involve all of the above behaviors. In fact, it is 
indeed possible to find quantum systems in which there is an abrupt change of the inner structure of 
the ground state that can be detected by some properties but not by others [17]. 

Let us give a simple example of a quantum critical point: consider the (1 + l)-dimensional a 
ferromagnetic quantum Ising spin chain, as defined by the Hamiltonian 



where cr a is the Pauli matrix a at site i of the chain, J > is a coupling parameter, and N is the 
number of spins. At J = oo the ground state of the system is two-fold degenerate and consists of 
all the spins aligned ferromagnetically in the jc-direction, being its subspace spanned by the two 
vectors |+, +, . . . , +) and |-, -, . . . , ->, where |+) and |-) denote the two possible eigenstates of the 
pauli matrix cr x . On the other hand, at J - the ground state of the system consists of all the spins 
aligned along the z-direction, |0, 0, . . . , 0>, where |0> = + I - ))- We now consider the behavior 

of the magnetization per particle of the ground state in the z-direction, as defined by the expected 

value M = — '■ . In the thermodynamic limit N — > oo this quantity tends to one when ./ — > 0, 
and tends to zero when J — > oo. A detailed analysis of this model in this limit shows that there is 
a specific point at which the magnetization per particle has a sudden change, as is represented in 
FigH This behavior implies that the model undergoes a second-order quantum phase transition at 
the critical point / = 7* = 1 in the thermodynamic limit. 

One may wonder what is the symmetry that we are breaking in this simple example of a quantum 
phase transition: it is the symmetry Z2 that the Hamiltonian from Eq|24]has at high values of the 
coupling parameter. In fact, this symmetry could even be further broken when J — > 00 if some 
extremely small magnetic field in the x-direction were present in our system, selecting one of the 
two possible ground states within this phase. In such a case, it is said that the symmetry of the 
Hamiltonian is spontaneously broken. 

A useful tool in the study of quantum critical systems is the renormalization group [18, 19], 
which describes the way in which a theory gets modified under scale transformations. Given some 
Hamiltonian depending on a set of parameters, the transformations of the renormalization group 
define a flow in the parameter space, and in particular the fixed points of those transformations 
correspond to theories which are invariant under changes of scale. Indeed, the essence of the renor- 
malization procedure is the elimination of degrees of freedom in the description of a system. This 
point of view is one of the basis for the development of different numerical techniques that allow 
to compute basic properties of quantum many-body systems, as is the case of the so-called density 
matrix renormalization group algorithm [20]. 

The behavior of many quantum critical models can also be explained by using tools from con- 
formal field theory [21]. There are quantum many -body systems which can be understood as a 
regularization on a lattice of a quantum field theory, as is the case of the previously-discussed Ising 



N 



N 




(24) 



1=1 



1=1 



a We use the field-theoretical notation (1 + 1) to denote one spatial and one temporal dimension. Time is always to be 
kept fixed. 
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Figure 3: Magnetization per particle in the ferromagnetic quantum Ising spin chain as a function 
of the coupling parameter, in the thermodynamic limit. The point J - J* - I corresponds to a 
second-order quantum phase transition point. 

quantum spin chain, which can be represented by the quantum field of a (1 + l)-dimensional spin- 
less fermion [22]. When those quantum many-body systems become critical, their description in 
terms of a quantum field theory allows to see that the symmetry group is not composed of only 
scale transformations, but of the full group of conformal transformations. In fact, conformal sym- 
metry is particularly powerful when applied to (1 + l)-dimensional quantum systems, allowing to 
determine almost all the basic properties of the model in consideration just by means of symmetry 
arguments. We perform some conformal field theory calculations in this thesis, and some basic 
technical background is given in Appendix B. 

Experimental quantum computers 

There will exist some day a quantum computer? This apparently simple question is by no means 
easy to answer. Actually, it is the opinion of some scientists that it is eventually impossible to build 
a quantum computer because of the unavoidable problem of the decoherence that any quantum 
system undergoes when it interacts with its environment. Nevertheless, other physicists think that 
these experimental drawbacks can be eventually in part ameliorated if the appropriate conditions 
are given. The main requirements that any experimental proposal must match if its purpose is to 
faithfully represent a quantum computer are known as the DiVincenzo criteria [23], and so far there 
have been many different ideas to perform experimental quantum computation that try to fulfill as 
much as possible these conditions. Important proposals are those based on quantum optical devices, 
such as the optical photon quantum computer, cavity quantum electrodynamics devices, optical 
lattices, or ion traps [24]. The idea of performing quantum computation by means of nuclear 
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magnetic resonance (NMR) has been considered as well [25-27]. Furthermore, proposals based 
on superconductor devices, quantum dots [28], and doped semiconductors [29,30] have also been 
considered by different people. The future development of these and other experimental techniques, 
and to what extent they can implement a many-qubit quantum computer, remains yet uncertain. A 
detailed discussion about experimental quantum computation can be found for instance in [2]. 

What is this thesis about? 

We focus here on the fields of quantum information science, condensed-matter physics, and quantum 
field theory. While these three branches of physics can be regarded as independent by themselves, 
there are clear overlaps among them, such that knowledge from one field benefits the others. As 
we said, conformal field theory [21] has helped to understand the universality classes of many crit- 
ical (1 + l)-dimensional quantum many-body systems. Also, the study of the entanglement present 
in the ground state of quantum Hamiltonians at a quantum phase transition shows direct analogies 
with those coming from the study of entropies in quantum field theory [31^14]. These results in 
turn connect with the performance of numerical techniques like the density matrix renormalization 
group [20], that allow to compute basic properties of some quantum many -body systems [45-60]. 
Indeed, quantum phase transitions are very much related to the model of adiabatic quantum compu- 
tation [16,61-71], which poses today challenges within the field of computational complexity [72]. 

The work that we present in this thesis tries to be at the crossover of quantum information 
science, quantum many -body physics, and quantum field theory. We use tools from these three fields 
to analyze problems that arise in the interdisciplinary intersection. More concretely, in Chapter 1 we 
consider the irreversibility of renormalization group flows from a quantum information perspective 
by using majorization theory and conformal field theory. In Chapter 2 we compute the entanglement 
of a single copy of a bipartite quantum system for a variety of models by using techniques from 
conformal field theory and Toeplitz matrices. The entanglement entropy of the so-called Lipkin- 
Meshkov-Glick model is computed in Chapter 3, showing analogies with that of (1 + l)-dimensional 
quantum systems. In Chapter 4 we apply the ideas of scaling of quantum correlations in quantum 
phase transitions to the study of quantum algorithms, focusing on Shor's factorization algorithm and 
quantum algorithms by adiabatic evolution solving an NP-complete and the searching problems. 
Also, in Chapter 5 we use techniques originally inspired by condensed-matter physics to develop 
classical simulations, using the so-called matrix product states, of an adiabatic quantum algorithm. 
Finally, in Chapter 6 we consider the behavior of some families of quantum algorithms from the 
perspective of majorization theory. 

The structure within each Chapter is such that the last section always summarizes the basic res- 
ults. Some general conclusions and possible future directions are briefly discussed in Chapter 7. 
Appendix A, Appendix B and Appendix C respectively deal with some basic notions on majoriza- 
tion theory, conformal field theory, and classical complexity theory. 



Chapter 1 

Majorization along parameter and 
renormalization group flows 

Is it possible to somehow relate physical theories that describe Nature at different scales? Say, 
given a theory describing Nature at high energies, we should demand that the effective low-energy 
behavior should be obtained by integrating out the high-energy degrees of freedom, thus getting a 
new theory correctly describing the low-energy sector of the original theory. This should be much 
in the same way as Maxwell's electromagnetism correctly describes the low-energy behavior of 
quantum electrodynamics. 

This non-perturbative approach to the fundamental theories governing Nature was essentially 
developed by Wilson and is the key ingredient of the so-called renormalization group [18, 19, 73]: 
effective low-energy theories can be obtained from high-energy theories by conveniently eliminating 
the high-energy degrees of freedom. To be more precise, the renormalization group is the mechan- 
ism that controls the modification of a physical theory through a change of scale. Renormalization 
group transformations then define a flow in the space of theories from high energies (ultraviolet 
theories) to low energies (infrared theories). Actually, it is possible to extend this idea, and the 
renormalization procedure can be more generically understood as the elimination of some given de- 
grees of freedom which we are not interested in because of some reason. The name "renormalization 
group" is used due to historical reasons, since the set of transformations does not constitute a formal 
group from a mathematical point of view. 

Since the single process of integrating out modes seems to apparently be an irreversible op- 
eration by itself, one is naturally led to ask whether renormalization group flows are themselves 
irreversible. This question is in fact equivalent to asking whether there is a fundamental obstruction 
to recover microscopic physics from macroscopic physics, or more generally, whether there is a net 
information loss along renormalization group trajectories. While some theories may exhibit limit 
cycles in these flows, the question is under which conditions irreversibility remains. Efforts in this 
direction were originally carried by Wallace and Zia [74], while a key theorem was later proven by 
Zamolodchikov [75] in the context of (l+l)-dimensional quantum field theories: for every unitary, 
renormalizable, Poincare invariant quantum field theory, there exists a universal c-function which 
decreases along renormalization group flows, while it is only stationary at (conformal) fixed points, 
where it reduces to the central charge c of the conformal theory. This result sets an arrow on renor- 
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malization group flows, since it implies that a given theory can be the infrared (IR) realization of 
another ultraviolet (UV) theory only if their central charges satisfy the inequality cjr < cjjv- 

The following question then arises: "under which conditions irreversibility of renormalization 
group flows holds in higher dimensions?". This has been addressed from different perspectives 
[76-94]. It is our purpose here to provide a new point of view about this problem based on the 
accumulated knowledge from the field of quantum information science, by focusing first on the 
case of (1 + 1) dimensions. 

An important application of quantum information to quantum many-body physics has been the 
use of majorization theory [11-14] in order to analyze the structure present in the ground state - also 
called vacuum - of some models along renormalization group flows. Following this idea, in [95] it 
was originally proposed that irreversibility along the flows may be rooted in properties concerning 
only the vacuum, without necessity of accessing the whole Hamiltonian of the system and its full 
tower of eigenstates. Such an irreversibility was casted into the idea of an entanglement loss along 
renormalization group flows, which proceeded in three constructive steps for (l + l)-dimensional 
quantum systems: first, due to the fact that the central charge of a (l + l)-dimensional conformal field 
theory is in fact a genuine measure of the bipartite entanglement present in the ground state of the 
system [36-44], there is a global loss of entanglement due to the c-theorem of Zamolodchikov [75]; 
second, given a splitting of the system into two contiguous pieces, there is a monotonic loss of 
entanglement due to the numerically observed monotonicity for the entanglement entropy between 
the two subsystems along the flow, decreasing when going away from the critical fixed - ultraviolet 
- point; third, this loss of entanglement is seen to be fine-grained, since it follows from a strict set 
of majorization ordering relations, numerically obeyed by the eigenvalues of the reduced density 
matrix of the subsystems. This last step motivated the authors of [95] to conjecture that there 
was a fine-grained entanglement loss along renormalization group flows rooted only in properties 
of the vacuum, at least for (l+l)-dimensional quantum systems. In fact, a similar fine-grained 
entanglement loss had already been numerically observed in [37, 38], for changes in the size of the 
bipartition described by the corresponding ground-state density operators, at conformally-invariant 
critical points. 

The aim of this Chapter is to analytically prove relations between conformal field theory, renor- 
malization group and entanglement. We develop, in the bipartite scenario, a detailed and analytical 
study of the majorization properties of the eigenvalue spectrum obtained from the reduced density 
matrices of the ground state for a variety of (l + l)-dimensional quantum models. Our approach 
is based on infinitesimal variations of the parameters defining the model - magnetic fields, aniso- 
tropics - or deformations in the size of the block L for one of the subsystems. We prove in these 
situations that there are strict majorization relations underlying the structure of the eigenvalues of 
the considered reduced density matrices or, in other words, that there is a fine-grained entanglement 
loss. The result of our study is presented in terms of two theorems. On the one hand, we are able 
to prove continuous majorization relations as a function of the parameters defining the model under 
study. Some of these flows in parameter space may indeed be understood as renormalization group 
flows for a particular class of integrable theories, like the Ising quantum spin chain. On the other 
hand, using the machinery of conformal field theory in the bulk we are able to prove exact continu- 
ous majorization relations in terms of deformations of the size of the block L that is considered. We 
also provide explicit analytical examples for models with a boundary based on previous work of 
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Peschel, Kaulke and Legeza [96-98]. 

1.1 Global, monotonous and fine-grained entanglement loss 

Consider the pure ground state |Q) of a given regularized physical system which depends on a 
particular set of parameters, and let us perform a repartition of the system into two pieces A and 
B. The density matrix for A, describing all the physical observables accessible to A, is given by 
Pa = tr^OOXni) - and analogously for B -. Here we will focus our discussion on the density matrix 
for the subsystem A, so we will drop the subindex A from our notation. Let us consider a change in 
one of the parameters on which the resultant density matrix depends, say, parameter "f, which can 
be an original parameter of the system, or be related to the size of the region A. To be precise, we 
perform a change in the parameter space from t\ to t 2 , with t% > t\. This involves a flow in the space 
of reduced density matrices from p(?i) to p(t 2 ), as represented in Fig ll.ll 



We wish to understand how this variation of the parameter alters the inner structure of the ground 
state and, in particular, how it modifies the entanglement between the two partys, A and B. Because 
we are considering entanglement at two different points t 2 and t\ , let us assume that the entanglement 
between A and B is larger at the point t\ than at the point t 2 , so we have an entanglement loss when 
going from t\ to t 2 . 

Our characterization of this entanglement loss will progress through three stages, refining at 
every step the underlying ordering of quantum correlations. These three stages will be respectively 
called global, monotonous and fine-grained entanglement loss. 



Global entanglement loss.- A possible way to quantify the loss of entanglement between A and 
B when going from t\ to t 2 is by means of the entanglement entropy S{p{t)) = -tr(p(?)log 2 p(0)- 
Since at t 2 the two partys are less entangled than at t\, we have that 



9(h) 




9(h) 



Figure 1.1: A flow in the space of density matrices, driven by parameter t. 



S(p(h))>S(p(t 2 )), 



(1.1) 



which is a global assessment between points t 2 and t\ . This is what we shall call global entanglement 
loss. 
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Monotonous entanglement loss.- A more refined condition of entanglement loss can be obtained 
by imposing the monotonicity of the derivative of the entanglement entropy when varying the para- 
meter "f. That is, the infinitesimal condition 

S(p(t))>S(p(t + di)) (1.2) 

implies a stronger condition on the structure of the ground state under deformations of the parameter 
along the flow in t. This monotonic behavior of the entanglement entropy is what we shall call 
monotonous entanglement loss. 



Fine-grained entanglement loss.- When monotonous entanglement loss holds, we can wonder 
whether the spectrum of p(t) becomes more and more ordered as we change the value of the para- 
meter. It is then plausible to ask if it is possible to make stronger claims than the inequalities given 
by Eq ll.ll and Eq ll.2l and unveil some richer structure. The finest notion of reordering when chan- 
ging the parameter is then given by the monotonic majorization of the eigenvalue distribution along 
the flow. If we call p(t) the vector corresponding to the probability distribution of the spectrum 
arising from the density operator p(t), then the infinitesimal condition 

f5{t)<frt + df) (1.3) 

along the flow in t reflects a strong ordering of the ground state along the flow. This is what we 
call fine-grained entanglement loss, because this condition involves a whole tower of inequalities 
to be simultaneously satisfied. This Chapter is devoted to this precise majorization condition in 
different circumstances when considering (1 + l)-dimensional quantum systems. For background 
on majorization, see Appendix A. 



1 .2 Maj orization along parameter flows in ( 1 + 1 ) -dimensional quantum 

systems 

Our aim in this section is to study strict continuous majorization relations along parameter flows, 
under the conditions of monotonicity of the eigenvalues of the reduced density matrix of the vacuum 
in parameter space. Some of these flows indeed coincide with renormalization group flows for some 
integrable theories, as is the case of the Ising quantum spin chain. 

Before entering into the main theorem of this section, let us perform a small calculation which 
will turn to be very useful: we want to compute the reduced density matrix for an interval of length 
L of the vacuum of a conformal field theory in (1 + 1) dimensions - see Appendix B for background 
on conformal field theory -. With this purpose, let Zi(q) = g~ c ^ 12 tr (g (Lo+ ^ o) ) denote the partition 
function of a subsystem of size L [21,36], where q - e 2mT , t = (/7r)/(ln (L/n)), q being an ultraviolet 
cut-off, and Lq and Lq the Oth Virasoro operators. Let b = c/12bea parameter that depends on 
the central charge and therefore on the universality class of the model. The unnormalized density 
matrix can then be written as q~ b q^ >+Lo \ since it can be understood as a propagator and (Lq + Lq) 
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is proportional to the generator of translations in time - which corresponds to dilatations in the 
conformal plane - [21]. Furthermore, we have that 

tr(g (Lo+Io) ) = 1 + mq ai + n 2 q ai +■■■ , (1.4) 

due to the fact that the operator (Lq + Lq) is diagonalized in terms of highest- weight states \h, h): 
(Lo + Lo)\h, h) = (h + h)\h, h), with h > and h > 0; the coefficients a\,a 2 , . . . > 0, a{ + aj V/ + j 
are related to the eigenvalues of (Lq+Lq), and n\, n 2 , . . . correspond to degeneracies. The normalized 
distinct eigenvalues of pi - YTl^' '^ (Lo+Lo) are then given by 



1 


(1 + n\q a i + n 2 q ai 


+ • 


■•) 


q ai 






(1 + mq a i + n 2 q ai 


+ • 


■•) 








(1 + niq a ^ + n 2 q ai 


+ • 


■•) 



(1.5) 



We are now in conditions of introducing the main result of this section, which can be casted into 
the following theorem: 



Theorem 1.1: Consider a (1 + l)-dimensional physical theory which depends on a set of real 
parameters g = (g\,g2, ■ ■ •)» such that 

• there is a non-trivial conformal point g*, for which the model is conformally invariant in the 
bulk, 

• the deformations from g* in parameter space in the positive direction of a given unity vector 
e preserve part of the conformal structure of the model, that is, the eigenvalues of the generic 
reduced density matrices of the vacuum p(g) are still of the form given by Eq \1.5\ with some 
parameter-dependent factors q(g), for values of the parameters g = g* + ae, and 

• the factor q(g) is a monotonic decreasing function along the direction ofe, that is, we demand 
that 



dqiM) 
da 

along the flow. 



Then, away from the conformal point there is continuous majorization of the eigenvalues of 
the reduced density matrices of the ground state along the flow in the parameters g in the positive 
direction ofe (see Fis \1.2\ . that is, 

p{g\)<P{g2), 

gi = g + ae, g2 = g + a e, a > a . 
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Figure 1.2: A possible flow in the space of parameters in the direction of e. 



Proof: Let us define the quantity Z(q) = (1 + n^q ai + n2q" 2 + •••), where it is assumed that 
q = q(g), for values of g along the flow in a. Notice that at conformal points Z(q(g*)) is not invariant 
under modular transformations, as opposed to the partition function Z(q(g*)). The behavior of the 
eigenvalues in terms of deformations with respect to the parameter a follows from 

m = mz±*jL< , (1 . 8) 

da q da 

and therefore 

«L=*LLko. (1 .9) 

da da \Z{q) ) 

Because Ai is always the largest eigenvalue Va, the first cumulant automatically satisfies continuous 
majorization along the considered flow. The variation of the other eigenvalues Ai (I > 1) with respect 
to a reads as follows: 



dAj_ _d tq 



«(/-!) 



da da \ Z(q) 

q a ^~ l t Z(q)-l\dq 



z<„ r <w, -^r)^' (U0) 

Let us concentrate on the behavior of the second eigenvalue A2. We observe that two different 
situations can happen: 



if 

Z(q) 



ai _^LLi)> , (1.11) 



g ( M) ~ Z{ t I 1 I > V/ > 2 , (1.12) 



then since a(;-i) > a\ VZ > 2, we have that 

Z(q) 

which in turn implies that 

dA t 

- i <0VZ>2. (1.13) 
da 



1.2. Majorization along parameter flows in (1 + l)-dimensional quantum systems 17 



From this we have that the second cumulant satisfies 



da 



da 



U>2 ) 



>0, 



(1.14) 



thus fulfilling majorization. The same conclusion extends easily in this case to all the remain- 
ing cumulants, and therefore majorization is satisfied by the whole probability distribution. 



if 



then 



Z(q)-l \ „ 
a\ = < , 

Z(q) 



dM >Q 
da 



and therefore 



da 



>0, 



(1.15) 



(1.16) 



(1.17) 



so the second cumulant satisfies majorization, but nothing can be said from the previous three 
equations about the remaining cumulants. 



Proceeding with this analysis for each one of the eigenvalues we see that, if these are monotonically 
decreasing functions of a then majorization is fulfilled for the particular cumulant under study, 
but since or (+1 > or ; - Vz we notice that once the first monotonically increasing eigenvalue is found, 
majorization is directly satisfied by the whole distribution of eigenvalues, therefore p{g\) < p{g2) if 
gi - g* + ae, g2 - g* + a'e, and a' > a, as claimed. □ 



An interesting application of Theorem 1 . 1 comes whenever a can be related to the scale of a 
renormalization group transformation. Then it can be understood as a proof of fine-grained entan- 
glement loss along a renormalization group flow for a particular set of integrable theories, namely, 
those theories which fulfill the hypothesis of our theorem. We stress that, while it would probably 
be possible to obtain results based on perturbation theory in the neighborhood of the conformal 
point for non-integrable theories, our theorem is based on the alternative approach of completely 
non-perturbative results under the assumption of integrability of the theory along the flow. This 
assumption is naturally fulfilled by many interesting models: we wish to illustrate this point with 
the analytical examples of similar situations for the Heisenberg and XY quantum spin chains with 
a boundary. At this point we wish to remark as well that, for those theories depending only on 
one parameter g, the monotonicity in the change of the parameter along a renormalization group 
flow between two fixed points is trivial, since between two zeros the /J-f unction /3 = - jjfj, I being 
the scale of the renormalization group transformation, can only be either positive or negative, thus 
implying the monotonicity of the parameter when flowing from one fixed point to the other. Notice 
that our claim, which is majorization of the reduced density matrices of the vacuum, is stronger. 
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A majorization lemma 

As a previous step in our derivations, let us state a useful lemma about majorization theory which 
we shall constantly use in the forthcoming sections. We refer the reader to Appendix A for math- 
ematical definitions and more background on majorization theory. The lemma reads as follows: 

Lemma 1.1 [95]: If p\ < pi and q\ < q*2, then {p\ ® q\) < (p2 ® q2)- This means that 
majorization is preserved under the direct product operation. 

Proof: If pi < p*2 and q\ < q*2 then p\ = D p p2 and q\ - D q q*2 where D p ,D q are both doubly 
stochastic matrices. Therefore (p\®q\) = (D p ®D q )(p2®q2), where (D p ®D q ) is a doubly stochastic 
matrix in the direct product space, and so (p\ ®q\)<(p2® qi)- □ 

1.2.1 Quantum Heisenberg spin chain with a boundary 

Consider the Hamiltonian of the Heisenberg quantum spin chain with a boundary 

oo 

= (1.18) 

!=1 

where A > 1 is the anisotropy parameter. This model is non-critical in the region defined by A > 1 
and critical at A = 1 . Notice that, since this is a uniparametric theory which can be mapped to a 
Gaussian free theory, any renormalization group transformation must be reflected in a change of the 
only existing parameter. Thus, the flow in A must necessarily coincide with a renormalization group 
flow. 

From the pure ground state of the system, we trace out the Af/2 contiguous spins i - 1,2,..., N/2, 
getting an infinite-dimensional density matrix pa in the limit N — > oo which describes half of the 
system, and such that it can be written as a thermal density matrix of free fermions [96-98]. Its 
eigenvalues are given by 

p A («o, «!,...,««,) - ^-e~ nk£k 

= Pa(«o)Pa(«i)---Pa(«oo) , (1-19) 

with Pa(«/0 = where Z\ - (1 + e~ £k ) is the partition function for the mode k, n^ - 0, 1, for 

k - 0, 1 , . . . , oo and with dispersion relation 

e k - 2k arcosh(A) . (1.20) 

The physical branch of the function arcosh(A) is defined for A > 1 and is a monotonic increasing 
function of A. On top, the whole partition function Za can be decomposed as an infinite direct 
product of the different free fermionic modes: 



oo 

Z A = Y](l+e-*) 

k=0 



(1.21) 
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From the last equations, it is not difficult to see that pa -< pa' if A < A'. Fixing the attention 
on a particular mode k, we evaluate the derivative of the largest probability for this mode, P k A = 
(1 + e~ £k )~ . This derivative is seen to be 



dP k A 2k 



dA (1 + e -e k) i Va^H" 



>0, (1.22) 



for k - 1,2, ... oo and for k = 0. It follows from this fact that all the modes independently major- 
ize their respective probability distributions as A increases, with the peculiarity that the Oth mode 
remains unchanged along the flow, since its probability distribution is always (j, l) T . The particular 
behavior of this mode is responsible for the appearance of the "cat" state that is the ground state 
for large values of A - notice that in that limit the model corresponds to the quantum Ising model 
without magnetic field -. These results, together with the Lemma 1.1, make this example obey 
majorization along the flow in the parameter, which can indeed be understood as a renormalization 
group flow because of the reasons mentioned at the beginning of the example. 

1.2.2 Quantum XY spin chain with a boundary 

Similar results to the one obtained for the Heisenberg model can be obtained for a different model. 
Let us consider the quantum AT-model with a boundary, as described by the Hamiltonian 



(<I^ <1 + <Iz^ l + ^), (1 . 23) 



where y can be regarded as the anisotropy parameter and A as the magnetic field. The phase diagram 
of this model is shown in Fig ll .31 where one can see that there exist different critical regions de- 
pending on the values of the parameters, corresponding to different universality classes [37-40,99]. 
Similarly to the previous example, this model is integrable and can be mapped to a Gaussian free 
theory with a mass parameter depending on a particular combination of both A and y once the kinetic 
term has been properly normalized (see [22]). A renormalization group flow can then be understood 
as a set of flows in the plane of A and y. 

Consider the ground state of Eq l 1.231 and trace out the contiguous spins i = 1, 2, . . . , N/2 in the 
limit N — > oo. The resulting density matrix P(a, 7 ) can be written as a thermal state of free fermions, 
and its eigenvalues are given by [96-98]: 

1 _v~ 

P(A,y)(no,n u . ..,««,)= e ^ n ^ k > (1 24 ) 

where = 0, 1 , and the single-mode energies are given by 

[2ke , if A < 1 

«* = \ (1.25) 
\{2k+\)e, iO>l, 

with k = 0, 1, . . . , oo. The parameter e is defined by the relation 

/(VI-jc 2 ) 

e = n y ' , (1.26) 

I(x) 
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XX Ising 
critical XY 



critical Ising 



critical XX 







1 Y 

Figure 1.3: Phase diagram of the quantum AT-model. 



I(x) being the complete elliptic integral of the first kind 

•* /2 dG 



/(x)= f , d.27i 
Jo 



^1 - x 2 sin 2 (G) 



and x being given by 



T J(V^ 2 + y 2 -i)/y , if^<i 

\y/(V^ 2 + r 2 - 1), iO>l, 

where the condition A 2 + y 2 > 1 is assumed for a correct behavior of the above expressions (external 
region of the Baruoch-McCoy circle [99]). 

We observe that the probability distribution defined by the eigenvalues of P(A, y ) is again the 
direct product of distributions for each one of the separate modes. Therefore, in order to study 
majorization we can focus separately on each one of these modes, in the same way as we already 
did in the previous example. We wish now to consider our analysis in terms of the flows with respect 
to the magnetic field A and with respect to the anisotropy y in a separate way. Other trajectories in 
the parameter space may induce different behaviors, and a trajectory-dependent analysis should then 
be considered for each particular case. 



Flow along the magnetic field A 

We consider in this subsection a fixed value of y while the value of A changes, always fulfilling 
the condition A 2 + y 2 > 1. Therefore, at this point we can drop y from our notation. We separate 
the analysis of majorization for the regions 1 < A < oo and + -^1 - y 2 < A < 1 for reasons that 
will become clearer during the study example but that already can be realized just by looking at the 
phase space structure in Fig |1.3l 
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Region 1 < A < oo.- We show that px < px & A < A' .\xv this region of parameter space, the largest 
probability for the mode k is P* = (1 + e -< *) -1 . The variation of P*j with respect to A is 

dPj = (2k + l)e~^ de 
dA (i + e -(2k+\)e^- dA ' 

^ ^ — — — — - dP^ 

A direct computation using Ea fl~26l Ea lHTl and EqE33]shows that ^ > 0. Therefore, -gf > for 
k - 0, 1, . . . , oo. This derivation shows mode-by-mode majorization when A increases. Combining 
this result with the Lemma 1 . 1 , we see that this example obeys majorization. 



Region + ^1 — y 2 < A < 1.- For this case, we show that px < Px 1 if /I > A'. In particular, the 
probability distribution for the 0th fermionic mode remains constant and equal to (^, , which 
brings again a "cat" state for low values of A. Similarly to the latter case, the largest probability for 
mode k is P k A - (1 + e~ £k )~ l , with 

e k = 2kn ^ X ~ xl - = 2ke , (1.30) 
I(x) 

and x = ( sjA 2 + y 2 — l)/y. Its derivative with respect to /I is 

d/ 5 ^ 2fe" 2fc de 



dA + e -2fe) 2 d/l 



(1.31) 



It is easy to see that this time gf < 0, and therefore -jj- < for k = 1, 2, . . . , oo, which brings 
majorization individually for each one of these modes when A decreases. The mode k = calls 

dP* =0 

for special attention. From Eq ll.311 it is seen that —£- = 0, therefore the probability distribution 
for this mode remains equal to (^, all along the flow. This is a marginal mode that brings the 
system to a "cat" state that appears as ground state of the system for low values of A. Notice that 
this peculiarity is rooted on the particular form of the dispersion relation given in Eq ll.251 which is 
proportional to 2k instead of 2k + 1 for this region in parameter space. These results, together with 
the Lemma 1.1, prove that this example also fulfills majorization. 

Flow along the anisotropy y 

In this subsection, the magnetic field A is fixed and the anisotropy y is the only free parameter of the 
model, still fulfilling A 2 + y 2 > 1. Thus, at this point we can drop A from our notation. We will see 
that p y < p y > if y > y', in the two regions 1 < A < oo and + yjl - y 2 < A < 1. In particular, in the 
region + yjl - y 2 < A < 1 , the probability distribution for the 0th fermionic mode remains constant 
and equal to (^, . Let us consider the biggest probability for the mode k, P y = (l+ e~ £k )~ l , with 
ek = we, where 

[2k, if A < 1 

oj = I (1.32) 
(2£ + l), ifA>\, 
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and e as defined in the preceding sections. It is easy to verify that 



dy 



(1 + e-««*)2 dx dy 



aje-^ k de dx 



<0 



(1.33) 



for k - 0, 1, . . . , oo if A > 1 and for k - 1, 2, . . . , oo if A < 1. The mode k - for A < 1 needs of 



special attention, since —4— = 0, and therefore the probability distribution for this mode remains 
constant and equal to (±, ±) r all along the flow. These results, together with the Lemma 1.1, show 
that this case obeys again majorization along the flow in the parameter. 

1.3 Majorization with L in ( 1 + 1) -dimensional conformal field theories 

A similar study to the one presented in the previous section about majorization along flows in 
parameter space can be now performed exclusively at the conformal point for flows in the size of 
the block under consideration. Here we present an analytical derivation of majorization relations for 
any (1 + l)-dimensional conformal field theory without boundaries - or in the bulk a - in the bipartite 
scenario when the size of the considered subsystems changes, that is to say, under deformations in 
the interval of the accessible region for one of the two partys. This size will be represented by the 
length L of the space interval for which we consider the reduced density matrix after tracing out 
all the degrees of freedom corresponding to the rest of the universe. Our main result in this section 
can be casted into the following theorem: 

Theorem 1.2: pi < pi> if L > L' for all possible (1 + l)-dimensional conformal field theories 
in the bulk. 

Proof: Since the factors q are now monotonic functions of the size of the interval L, the proof 
of this theorem is analogous to the proof of Theorem 1.1, with the only exception that now the cu- 
mulants are monotonically decreasing (instead of increasing) functions along the flow in L. Taking 
this into account, it immediately follows that pi < pu if L > U . This proof is valid for all possible 
(1 + l)-dimensional conformal field theories in the bulk, since it only relies on completely general 
assumptions. □ 

1.3.1 Critical quantum XX spin chain with a boundary 

Let us give an example of a similar situation to the one presented in Theorem 1 .2 for the particular 
case of the quantum XX-mode\ with a boundary, for which the exact spectrum of pr, can be explicitly 
computed. The Hamiltonian of the model without magnetic field is given by 



dP*=° 



OO 




(1.34) 



i=\ 



a The case in which boundaries are present in the system must be properly considered from the point of view of the 
so-called boundary conformal field theory. This has been done by H.Q. Zhou et al. in [100]. For technical background 
on conformal field theory without boundaries, see Appendix B. 
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The system as described by this model is critical since it is gapless. Notice that the ultraviolet cut-off 
coincides with the lattice spacing and the theory is naturally regularized, hence rj = 1 . Taking the 
ground state and tracing out all but a block of 1, 2, . . . ,L contiguous spins, the density matrix pi 
describing this block can be written, in the large-L limit, as a thermal state of free fermions [96-98]: 



e~ H ' 



PL = (1-35) 

Zl 



Zl being the partition function for a given L, H' = 11^=0 e kd[dk, with fermionic creation and anni- 
hilation operators ct k , d\ and dispersion relation 

n 2 

e k = ——(2k + l) k = 0,l,...,L-l. (1.36) 
2 InL 

The eigenvalues of the density matrix pi can then be written in terms of non-interactive fermionic 
modes 

p L (n ,n u ...,n L -i) = i-g" 

Zl (1-37) 

= pdm) • • • pdnL-i) , 

with p(nk) - ^e~ nk£k , where z£ = (1 + e~ £k ) is the partition function for the mode k, and nt - 0, 1, 
Vk. It is worth noticing that the partition function of the whole block Zl factorizes as a product over 
the L modes: 

L-l 

Z L = ]~[(l+0 ■ (1.38) 

k=0 

Once the density matrix of the subsystem is well characterized with respect to its size L, it is 
not difficult to prove that pl < pi> if L > L' . In order to see this, we will fix our attention on 
the majorization within each mode and then we will apply Lemma 1.1 for the whole subsystem. 
We initially have to observe the behavior in L of the largest probability defined by each individual 
distribution for each one of the modes, that is, P k L - 1 /z£ = (1 + e~ ek )~ l , for k - 0, 1, . . . , L - 1. It 
is straightforward to see that 



AL (i + e -tkf dL 



<0, (1.39) 



which implies that P k L decreases if L increases This involves majorization within each mode 
k - 0, 1 , . . . , L - 2 when decreasing L by one unit. In addition, we need to see what happens with 
the last mode k - L - 1 when the size of the system is reduced from L to L - 1. Because this 
mode disappears for the system of size L - 1, its probability distribution turns out to be represented 
by the probability vector (l,0) r , which majorizes any probability distribution of two components. 
Combining these results with Lemma 1.1, we see that this example for the quantum XX-model 
provides a similar situation for a model with a boundary to the one presented in Theorem 1.2. 
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1.4 Conclusions of Chapter 1 

In this Chapter we have analyzed majorization relations along parameter and renormalization group 
flows for a variety of models in (1 + 1) dimensions. We have also provided in a rigorous way 
explicit and detailed proofs for all the majorization conjectures raised in some papers on quantum 
spin chains [37,38,95]. In order to be more specific: 

• We have proven the existence of a fine-grained entanglement loss for (1 + l)-dimensional 
quantum systems along uniparametric flows, when perturbations in parameter space preserve 
part of the conformal structure of the partition function, and some monotonicity conditions 
hold as well. These flows may coincide with renormalization group flows in some cases. We 
also considered similar situations which can be treated analytically, arising in the Heisenberg 
and XY models with a boundary. 

• We have also developed a completely general proof of majorization relations underlying the 
structure of the vacuum with respect to the size of the block L for all possible (1 + 1)- 
dimensional conformal field theories in the bulk. An example of a similar situation has been 
considered for the particular case of the ZZ-model with a boundary. 

These results provide solid mathematical grounds for the existence of majorization relations 
along renormalization group flows underlying the structure of the vacuum of (1 + l)-dimensional 
quantum spin chains. It would be interesting to relate the results of this Chapter to possible exten- 
sions of the c-theorem [75] to systems with more than (1 + 1) dimensions. While other approaches 
are also possible [76-88], majorization may be a unique tool in order to assess irreversibility of 
renormalization group flows in terms of properties of the vacuum only, and some numerical results 
in this direction have already been observed in systems of different dimensionality for flows in the 
parameter space [101, 102]. The analytical derivation and the consideration of the consequences for 
higher-dimensional systems of the properties presented here for (1 + 1) dimensions remains an open 
problem. 



Chapter 2 



Single-copy entanglement in 

(1 + 1) -dimensional quantum systems 

How much entanglement is contained in a given quantum many-body system? This simple but fun- 
damental question has been considered for systems close to and at quantum phase transitions by 
means of analyzing very different entanglement measures [17,31,36-44,56, 103-118]. All these 
different ways of measuring entanglement lead to results which complement each other and which 
help us to understand the precise way in which the ground state of critical models is organized. 
While the concurrence measures the pairwise entanglement that is present in the system between 
two of its specific constituents [119], the entanglement entropy measures the entanglement that ap- 
pears between two different blocks in a bipartition, in turn showing very interesting connections to 
the entropic area law found for systems such as black holes [31-35]. A detailed analysis of the en- 
tanglement entropy in critical quantum spin chains unveils a universal logarithmic scaling law with 
the size of the block under consideration, which admits an explanation in terms of the underlying 
conformal field theory in (1 + 1) dimensions [36-44]. Furthermore, it is now well understood that 
the good performance of density matrix renormalization group algorithms in (1 + 1) dimensions 
relies very much on this property 3 [56]. 

Our aim in this Chapter is to study an entanglement measure which, very much like the entan- 
glement entropy, is proven to have intriguing scaling properties for (1 + l)-dimensional quantum 
systems. We call this measure single-copy entanglement [113, 120], and its operational definition 
comes naturally motivated by a practical reason: while the entanglement entropy measures the aver- 
age amount of entanglement possible to be distilled from a bipartite system in the limit of having an 
infinite number of copies of the system [121], the single-copy entanglement measures the amount 
of entanglement present in the more realistic case of having just one copy of the system, in a way 
to be precisely defined later. As we shall see, we are able to analytically compute the asymptotic 
leading scaling behavior of the single-copy entanglement for all (1 + l)-dimensional conformal field 
theories in the bulk, together with its first-order correction. At that point in our derivations a sur- 
prise will appear: the entanglement contained in a single specimen of a critical (1 + l)-dimensional 
system is seen to be, asymptotically, half the entanglement that is available in the ideal case of 

"The relation between scaling of entanglement and the performance of classical numerical simulations for different 
quantum systems will be addressed in detail in Chapters 4 and 5. 
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having an infinite number of copies. This result is reinforced by an analysis from the point of view 
of quasi-free fermionic systems in (1 + 1) dimensions which leads again to similar conclusions: 
whenever the entanglement entropy scales logarithmically in the size of the system, the single-copy 
entanglement scales asymptotically as half of the entanglement entropy. Furthermore, and in order 
to make our study more complete, we also analyze the behavior of single-copy entanglement away 
from criticality for the specific example of the XY quantum spin chain. Let us then begin our study 
by formally defining what the single-copy entanglement is. 



2.1 Operational definition of the single-copy entanglement 

Let us ask ourselves the following question: how much entanglement is contained in an infin- 
ite number of copies of a pure bipartite system Ii^ab)? Let us be more specific with the term 
"how much", by posing the question differently: what is the maximal rate at which EPR-pairs 
-4= (|0)a|0)b + |1>a|1)b) can be distilled from an infinite number of copies of a pure bipartite sys- 
tem |(Aab)> just by invoking local operations and classical communication (LOCC) between the two 
parties? The answer to this question was originally found by Bennett et al. in [121]: if we are able 
to distill M EPR-pairs from N copies of a pure bipartite system |(Aab)> the rate M/N coincides, in 
the infinite-copy limit, with the entanglement entropy between the two partys, namely 

M 

hm — = S(pa) = -tr(p A log 2 p A ) = S(p B ) = -tr(p B log 2 p B ) , (2.1) 

JV->oo Jy 

Pa and p# respectively being the reduced density matrices of the two partys A (Alice) and B (Bob). 
This situation corresponds to the one represented in Fig l2.ll 

Alice QoD 





Figure 2.1: Scenario defining the entanglement entropy. Alice and Bob share an infinite number of 
copies of the bipartite system, and wish to distill EPR-pairs by performing LOCC. 
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While the above definition of entanglement entropy obviously makes sense, having an infinite 
number of copies of the system at hand is an unrealistic situation from the experimental point of 
view. Thus, let us now ask ourselves this variant of the above original question: how much entan- 
glement is contained in a single specimen of a pure bipartite system |^ab)? Or, equivalently, what 
is the largest entanglement content that any apparatus could potentially distill by LOCC from just 
one bipartite entangled system at hand? This scenario is represented in Fie l2.2l 

Alice Bob 





Figure 2.2: Scenario defining the single-copy entanglement. Alice and Bob share only one copy of 
the bipartite system, and wish to distill a maximally entangled state of the largest possible dimension 
by performing LOCC. 

The maximum entanglement that it is possible to obtain by distillation with LOCC in the single- 
copy case can be measured by the largest dimension of a maximally entangled state that can be 
distilled with certainty from the single specimen. That is, for a pure bipartite state Ii^as) with 
reduced density matrices p A and p#, we write for the single-copy entanglement 

Ei (p A ) = £i(p B ) - log 2 (M) (2.2) 

if 

I<Aab> i— > I<Am> under LOCC , (2.3) 

where 

1 M 

|<Am> = -= J] \q) A \q)B (2.4) 

is a maximally entangled state of dimension M. Now, we recall the result that the interconversion of 
bipartite pure states under LOCC in the single-copy case is governed by the following majorization 
relation for the reduced density matrices [15]: 

I<Aas> — > I«Aab> under LOCC <^> p A <p A , (2.5) 
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where pa is the reduced density matrix of the converted state \4>ab) for the party A b . Replacing in the 
above condition \i//ab) - I<Am) an d Pa = m^-m, Im being the M X M identity matrix, and considering 
the definition of majorization between probability distributions in terms of a set of inequalities to be 
satisfied by partial sums of its components - see Appendix A -, we find the inequality 

h < -J- => M < -J- , (2.6) 

Ai being the largest eigenvalue of pa- Given the above upper bound for M, one finds that 

£i (p A ) = -log 2 ii =Ei(p B ). (2.7) 

Therefore, the single-copy entanglement can be directly computed by looking only at the largest 
eigenvalue of the reduced density matrix of the system under consideration. This situation is very 
different from that of the entanglement entropy, where all the eigenvalues of the reduced density 
matrix contribute to the final quantity. 



2.2 Exact conformal field theoretical computation 

Now we wish to show the exact and analytical computation of the single-copy entanglement in 
the case of (1 + l)-dimensional conformal field theories in the bulk. We remind that the systems 
described by these theories correspond to the continuum limit of a variety of regularized quantum 
critical theories defined on a chain. For technical background, see Appendix B. 

As we saw in the previous Chapter, the reduced density matrix for a block of size L describing 
the vacuum of a (1 + l)-dimensional conformal field theory can be written as [21,36, 112] 

PL = -L- q -^2 q (L 0+ L 0) ^ (2 g) 

Zdq) 

where c is the central charge of the theory, Lq and Lq are the Oth holomorphic and antiholomorphic 
Virasoro operators, Zdq) is the partition function, q - e 2mT , and r = (in)/(\n(L/r])), rj being a 
regularization ultraviolet cut-off. For critical quantum chains we have that r] - 1 , which corresponds 
to the lattice spacing, and which is to be understood in our forthcoming calculations. 

The largest eigenvalue of the density matrix pi corresponds to the zero mode of (Lo + lo), that 

is, 

* = T7-^ ClU ' (2 " 9) 

since for this mode |0) we have that (Lo + Lo)|0) = 0. We then get a first expression for the single- 
copy entanglement: 

E 1 (p L ) = -\og 2 A l =\og 2 (Zdq)q cn2 ). (2.10) 

The leading behavior for the partition function can be computed when L is large by taking advantage 
of its invariance under modular transformations. The needed transformation corresponds to t — > 
-1/t, which amounts to Zdq) = Zdq), q = «r 27r2 / lnL , q = e~ 2lnL = 2~ 21o ^ L . It is now possible to 



b Of course the same relation holds as well for the party B. 
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expand the partition function in powers of q, since all the eigenvalues of the operator (Lq + Lq) are 
positive, and find that the leading contribution originates from the central charge: 

\og 2 Zdq) = ~\og 2 q + o( l -\= C -\og 2 L + o(\\ . (2.11) 



12 

This result translates into an explicit expression for the single-copy entanglement 

^ = l^-i^r +0 (z)- 

We wish to point out that the above result is exact up to polynomial corrections in 1 /L since no 
further powers of 1 / log 2 L appear in the expansion when L is large. 

Similar conformal field theory manipulations were used to prove that the von Neumann entropy 
for the same reduced density matrix pi is given by [36] 

S(p L ) = -^log 2 <? + Ofij , (2.13) 

which implies the following direct relation between entropy and single-copy entanglement: 

1 c (n log 9 e) 2 I log 9 L \ 

2 6 log 2 L \ L ) 

where the last subleading correction is easily calculated by comparing the results from [36] and 
our expression given in Eq l2.12l It should be noted here that the above result completely fixes the 
leading eigenvalue of the reduced density matrix of the block of size L to be dictated by its entropy 
within the large-L limit, that is, 

M 



lim — — = 1 . (2.15) 

Collections to this limit can be obtained from Eq l2.14l Quite remarkably, we also notice that all 
the eigenvalues will inherit the same leading behavior and differ by their subleading corrections 
controlled by the conformal weights corresponding to the universality class of the particular model 
in consideration. 



2.3 Exact computation in quasi-free fermionic quantum spin chains 

We aim now to reinforce the previously achieved result by investigating the same question from an 
alternative point of view, namely, we investigate all translationally invariant quantum spin models 
which can, under a Jordan-Wigner transformation, be written as an isotropic quadratic Hamiltonian 
in fermionic operators. 

The Jordan-Wigner transformation relates the Pauli operators in the quantum spin system to 
spinless fermionic operators {c ; } obeying the fermionic anticommutation relations 

{cj,c k } = 
[c),cl) = 

{c),c k \ = 6j k , (2.16) 
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according to 

1 /_1 

1 n=l 
1 M 

n=l 

a? - (2.17) 

Consider now an infinite quantum spin system in (1 + 1) dimensions that corresponds to a general 
translationally invariant isotropic quasi-free fermionic model. These correspond to chain systems 
whose Hamiltonian can be cast into the form 

H = Y,c]A M (2.18) 

l,k 

with A{ = A_i e R. The ground state of H is a quasi-free fermionic state, that is, a state that is 
completely characterized by the second moments of the fermionic operators. Notice that, while 
some of the spin chains described by this setting can be considered as well within the framework of 
conformal field theory in (1 + 1) dimensions, there may also be models that do not correspond to 
any such conformal field theory. 

Our claim is the following: if the entropy of entanglement satisfies 



S(p L )=flog 2 (L) + 0(l), (2.19) 
for some £ > 0, then the single-copy entanglement satisfies 

Ei(p L ) = \s(p L ) + 0(\). (2.20) 

That is, if we find that the entropy of entanglement scales asymptotically as the logarithm of L - 
as typically observed for this class of systems at criticality - then we can infer that the leading 
behavior of the single-copy entanglement will asymptotically be exactly one half of it. Notice that 
this does not fix such a relationship in the case that, for example, the system is gapped and the 
entropy of entanglement saturates (we shall consider an example of non-critical behavior within the 
next section). Let us now show how we arrive to the previous statement. 

The reduced state of a block of length L is entirely specified by the eigenvalues of the real 
symmetric Lx L Toeplitz matrix Tl, with Z-th row being given by (f_; + i, ?_/+2, to, tL-i)- The 
latter numbers are for an infinite quasi-free fermionic quantum chain found to be 

ti = ^ f 2 " g(k)e- m dk , (2.21) 

where g : C — > C is the so-called symbol [99, 122, 123], which essentially characterizes the fermi- 
onic model. The fact that Tl is a Toeplitz matrix reflects the translational invariance of the model. 
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The real eigenvalues of Ti will be labeled as m, ...jii e [-1, 1]. They can be found from the zeroes 
of the characteristic polynomial F : C — » C, 



F(z) = det(zl L -T L ) . 
The entropy of entanglement can then be obtained as [39,40, 108, 115] 



S(Pl) = J/s(1,w), 



where f s : R + x C —> C is defined as 

/*(*,?) = -(^) log; 

In fact, we can write [39,40, 108, 115] 



/=i 



x + y\ l x ~y\-, i x ~y 

'lor 



F'{z) 



S (p L ) = lim Um f f s (l+e, z) 1 =^dz 
s->0 s->0 2m J F(z) 



(2.22) 



(2.23) 



(2.24) 



(2.25) 



The contour of the integration in the complex plane is shown in Fig l2.3l In turn, we may also write 
for the single-copy entanglement [113] 



£i(Pl) = 2/i(0,ai/), 



(2.26) 



/=i 

in terms of the above fi\, where now f\ : R + x C — > C is to be defined as 

, , / 1 + (Z 2 +£ 2 ) 1/ 2 ' 
fi(s,z) - -log 2 l 

Respecting the cuts of the logarithm (see [113]), we may now cast E\(p£) into the form 

id 



(2.27) 




-id 




Figure 2.3: Contour of integration to be taken in case of both the entropy of entanglement and the 
single-copy entanglement. 



Ei(Pl) = lim lim -L [ f l{s ,z)^-dz . (2.28) 

Now we take advantage of the fact that Tl is a real symmetric Toeplitz matrix, which means that 
we can assess the asymptotic behavior of their determinants using proven instances of the Fisher- 
Hartwig conjecture [39,40,99, 108, 115, 122, 123]. We wish to remark at this point that the obser- 
vation that we only refer to proven instances of the Fisher-Hartwig conjecture derives from the fact 
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that we are only considering isotropic models [108]. Concerning the function F : C — » C, the 
Fisher-Hartwig conjecture allows us to write 

- a{z)L - b(z) log 2 L + 0(1) (2.29) 

in the large-L limit, where 

R 

b(z) = -2^(z)P'(z), (2-30) 

r=l 

with fi : C -» C such that [108] 

The number /?, in turn conesponds to half the number of discontinuities of the above symbol g(k) in 
the interval [0, 2n). Now, if we assume the validity of the logarithmic scaling of the entropy given 
in the expression of Eq l2.19l we know that, necessarily, 



lim lim f f s (l + s,z)a(z)dz = , (2.32) 

£->0 <5->0 J 

since no linear dependence in L must appear. Moreover, we know that S(pl) > E\(p£), which can 
easily be proven from their respective mathematical definitions - apart from the intuition that many 
copies of a system may help in entanglement distillation -. Therefore, in the large-L limit we must 
also necessarily have 

lim lim | Ms, z)a(z)dz - . (2.33) 

e^O 6->0 J 

Consequently, we only have to consider the logarithmically divergent term. For the entropy of 
entanglement the only relevant contour integral reads 

Is = lim lim — f f s (l + s,z)b(z)dz . (2.34) 
£->0 c5->0 2ni J 

In turn, for the single-copy entanglement the relevant contour integral becomes 

h = lim lim -L f Ms,z)b(z)dz . (2.35) 
£^0 6-*0 2ni J 

Taking into account that b(z) is analytic outside the interval [—1,1], the contributions of the circle 
pieces vanish in the two cases. Hence, we finally arrive at 

Sip L ) = 4 f dx^^log 2 (L) + 0{l) 
n L J_i 1 - x l 

R f 1 fi(0,x) 

Ei(p L ) = rfx^ylog 2 (L) + 0(l). (2.36) 
n L J_i 1 - x l 

Since /i(0, x) — — log 2 ((l + \x\)/2) for x e [-1, 1], we have that within the large-L limit, 

SiPL) = |log 2 L + 0(l) 

Ei(p L ) = ^log 2 L + 0(l), (2.37) 
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which in turn implies the validity of the expression that we anticipated in Eq 12.201 We have therefore 
proven that, in this class of models, whenever the system has a logarithmic asymptotical scaling of 
the entanglement entropy, the single-copy entanglement is exactly half the asymptotically available 
in the infinite-copy case in its leading contribution. We wish to remark as well that, from Eg 12.371 
the number R precisely corresponds to the central charge c for those models that are governed by 
an underlying conformal symmetry. For instance, for the quantum XX spin chain, we have that 
R = c - 1 , corresponding to the universality class of a free boson. 



2.4 Single-copy entanglement away from criticality 

In this section we exhibit an explicit example for which the relation between single-copy entangle- 
ment and entanglement entropy can be demonstrated near but off the critical region. We consider 
the XY quantum spin chain with a boundary, with Hamiltonian 

studied in Chapter 1 . Again, we consider the chain of semi-infinite length with a boundary, where 
the spins i - 1,2,..., N 72 with N — > oo have been traced out from the ground state of the system. 
The resultant density matrix P(x,y) can De written as a thermal density operator of a system of spinless 
fermions with creation and annihilation operators dt and d k in the following way [98]: 



e- H 



where 



H = Yt e kd\d k , (2.39) 



\2ke, ifA<\ 
e k = \ 2.40 
\(2k + 1)6 , if A > 1 , 

teN, and A e ft is the parameter controlling the external magnetic field, A* — 1 corresponding to 
the quantum phase transition point. We also have that 



„ /(VrT?) , (241) 
Ux) 



I(x) being the complete elliptic integral of the first kind, 

^ /2 de_ 

(l - x 2 sin 2 (6>)) 

Furthermore, x is related to the parameters A and y defining the model as follows: 



P (in 

Jo ^ • 1 - ] - 



f(V/l 2 + y 2 - f)/y . if/t<l, 
[y/(V^ 2 + 7 2 -l), iO>l, 



(2.43) 
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with the condition A 2 + y 2 > 1 (external region of the Baruoch-McCoy circle [99]). A computation 
of the single-copy entanglement with respect to this partitioning can be performed in terms of e, 
transforming sums into integrals by means of the Euler-McLaurin expansion, and finding 

7r 2 log 9 e e log 9 e 

EiifL^e) = — + 0{e~") (2.44) 

246 24 



if A < 1 and 



n log 9 e 1 e log 9 e 
Ei(p L ^ £ ) = + - + — jj- + 0(0 (2.45) 



if /i > 1 . No subleading corrections in powers of e do appear in the expansion. On the other hand 
it is easy to see by explicit evaluation that the entropy of entanglement can be related to the single 
copy-entanglement by 

S (p L ^oo,e) = (l-e|-Wl (p L ^oo,e) , (2.46) 



de 



which shows that 



lim [EipL^) - (pl^co,,)) - . (2.47) 



We notice that the limit e — > is precisely the limit where the theory becomes critical, that is when 
A — > A* - 1 . The above expression for finite e gives us corrections away from criticality to the 1 /2 
factor between the entanglement entropy and the single copy entanglement that has been discussed 
in the preceding sections. These corrections vanish as the system approaches criticality, as we have 
explicitly seen in this example. 



2.5 Conclusions of Chapter 2 

In this Chapter we have analyzed the single-copy entanglement, that is, the entanglement that it 
is possible to deterministically distill by using local operations and classical communication when 
only one copy of a bipartite system is at hand, in quantum systems in (1 + 1) dimensions. We have 
carried our analysis mainly from the point of view of conformal field theory in (1 + 1) dimensions 
in the bulk and quasi-free fermionic models in order to analyze critical systems, and also studied 
the behavior close to but away from criticality for the integrable example of the XY quantum spin 
chain. To be more precise: 

• For (1 + l)-dimensional conformal field theories we have proven that the leading scaling 
behavior of the single-copy entanglement is exactly half the asymptotic behavior of the en- 
tanglement entropy. The first-order correction to the leading term has also been explicitly 
computed. 

• For quasi-free fermionic quantum systems we have proven that if the asymptotic scaling of 
the entanglement entropy is logarithmic, then the asymptotic scaling of the single-copy entan- 
glement is also logarithmic, with a prefactor that is exactly half the one of the entanglement 
entropy. 
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• For the example of the semi-infinite XY quantum spin chain, we have computed the single- 
copy entanglement away from criticality and have observed that the factor 1 /2 between the 
entropy and the single-copy entanglement is only recovered when the system approaches the 
quantum phase transition point. 

The main conclusion is, therefore, that for (1 + l)-dimensional quantum systems at criticality 
the single-copy entanglement and the entanglement entropy for a system described by a reduced 
density matrix pi typically obey the law 

fm(|^ = 2). (2.48) 

For systems obeying the above relation we can say that in a single run, with a single invocation of 
a physical device acting on only one physical system, it is possible to obtain half the entanglement 
per specimen that is asymptotically available in the infinite-copy limit. Furthermore, all these res- 
ults also show relationships between the largest eigenvalue of the reduced vacuum p^ and its full 
spectrum for a very large class of quantum systems. 



Chapter 3 



Entanglement entropy in the 
Lipkin-Meshkov-Glick model 

Most of the analytical studies of the entanglement properties of quantum many-body systems close 
to criticality have been focused on the particular case of (1 + l)-dimensional systems, like the ones 
that we considered in the previous Chapters. Few models have been discussed so far in higher 
dimensions [31, 34, 35, 58-60, 101, 109-111, 124-131] either due to the absence of an exact di- 
agonalization of the system or to a difficult numerical treatment. We can in part understand this 
difficulty because of the existing link between the connectivity of a system and its entanglement 
entropy: one should naively think that, the bigger the connectivity of the system is, the bigger the 
amount of quantum correlations present in the ground state of the model should be, especially when 
the system is close to a quantum critical point. A classical numerical treatment of the model can 
become then very inefficient, as we shall in detail explain in the forthcoming Chapters 4 and 5. The 
idea in favor of this is rather simple: the more connected a system is, the more interactions it has, 
therefore the more entangled its ground state should be and the more difficult it should be to get its 
fundamental properties - like the ground-state energy or the correlation functions - by means of a 
classical numerical treatment. 

Actually, with some insight it is possible to make a non-accurate quantitative statement about 
the previous idea: given a system of N particles in (d + 1) dimensions, d being the number of spatial 
dimensions of the underlying lattice, if we believe that at criticality the entropy of entanglement 5 
is to scale proportionally to the area of the boundary of the region that separates the two subsystems 
under consideration, as is the case of bosonic systems [31, 105, 131], then it is not difficult to check 
that the entropy of a bipartition of the system between N/2 contiguous particles and the rest has to 
roughly scale like 

S-N^. (3.1) 

Critical fermionic systems may differ from the above law by means of an 0(log 2 AO multiplicative 
factor [109-111]. From the above reasoning we can see that the bigger the dimensionality d is - 
which is directly related to the connectivity of the system -, the stronger the scaling of the entan- 
glement entropy should be. The case of a conformally-invariant critical system with d = 1 has to 
be treated separately since the entropy has a logarithmic divergence, as we already remarked in pre- 
vious Chapters. This intuitive relation between entanglement and connectivity will be considered 
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again in Chapter 4, when studying the scaling of entanglement in quantum algorithms. 

In this context, the Lipkin-Meshkov-Glick model [132-134] has drawn much attention since 
it allows for a very efficient numerical treatment as well as analytical calculations. Furthermore, 
it provides a useful counter-example of the previous intuitive relation between entanglement and 
connectivity: in a system defined on a simplex - totally connected network -, and contrary to the 
intuition that we have specified before, the entanglement in the system behaves as if the system 
were (1 + l)-dimensional. This is a consequence of the role played by the symmetries within the 
description of the model, as we shall see. Entanglement can be increased by the connectivity, but 
can also be "killed" by the symmetries in some cases. 

First introduced by Lipkin, Meshkov and Glick in nuclear physics, this model has been the 
subject of intensive studies during the last two decades. It is of interest in order to describe in 
particular the Josephson effect in two-mode Bose-Einstein condensates [135, 136]. Its entanglement 
properties have been already discussed through the concurrence, which exhibits a cusp-like behavior 
at the critical point [137-139] as well as interesting dynamical properties [140]. Similar results have 
also been obtained in the Dicke model [141-143] which can be mapped onto the Lipkin-Meshkov- 
Glick model in some cases [144], or in the reduced BCS model [145]. Let us mention as well 
that the entanglement entropy has also been calculated for the anti-ferromagnetic Lipkin-Meshkov- 
Glick model [146] for which the ground state is known exactly [138, 147]. Here we analyze the 
von Neumann entropy computed from the ground state of the Lipkin-Meshkov-Glick model. We 
show that, at criticality, it behaves logarithmically with the size of the blocks L used in the bipartite 
decomposition of the density matrix with a prefactor that depends on the anisotropy parameter 
tuning the underlying universality class. We also discuss the dependence of the entropy with the 
magnetic field and stress the close analogy of the found results with those of (1 + l)-dimensional 
quantum systems. 



where cr? is the Pauli matrix at position k in the direction a, and N the total number of spins. This 
Hamiltonian describes a set of spins one-half located at the vertices of a Af-dimensional simplex - 
complete graph, as shown in Fig l3.ll - interacting via a ferromagnetic coupling A > in the xy- 
spin plane, y being an anisotropy parameter and h an external magnetic field applied along the z 
direction. 

Given that the model is defined on a simplex, the symmetry under permutations of particles 
allows us to rewrite the Hamiltonian from Eq l3.2l in terms of the total spin operators J a - ^i=i °^ /2. 
The previous Hamiltonian can then be expressed as 



3.1 The Lipkin-Meshkov-Glick model 



The Lipkin-Meshkov-Glick model is defined by the Hamiltonian 




(3.2) 



H = --(1 +y)(j 2 - J Z J Z -N/2)-2hJ z 



(3.3) 
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Figure 3.1: Complete graph - or simplex - of 8 vertices. 

where J 2 is the representation of spin N/2 of the Casimir operator and /* = J x ± iJ y . In the 
following, we set for simplicity A - 1 and since the spectrum of H is even under the transformation 
h <-> —h [140], we restrict our analysis to the region h > 0. Furthermore, we only consider the 
maximum spin sector J = N/2 to which the full spectrum of the Hamiltonian from Eq l4.18l belongs. 
A convenient basis of this subspace is spanned by the so-called Dicke states \N/2, M) which are 
invariant under the permutation of spins and are eigenstates of J 2 and J z with eigenvalues N(N+2)/4 
and M = -N/2, -N/2 + l,...,N/2- \,N/2, respectively. 



3.2 Entanglement within different regimes 

We consider the von Neumann entropy associated with the ground state reduced density matrix pl,n 
of a block of size L out of the total N spins, S l,n = S(pl,n) - -tr (plm^°&2Pl,n) and analyze 
its behavior as L is changed, both keeping N finite or sending it to infinity. Notice that since the 
ground state reduced density matrix is spanned by the set of (L + 1) Dicke states, the entropy of 
entanglement obeys the constraint S l,n ^ l°g2(^ + 1) f° r au ^ an d N, where the upper bound 
corresponds to the entropy of the maximally mixed state pl,n = I/(L + 1) in the Dicke basis. This 
argument implies that entanglement, as measured by the von Neumann entropy, cannot grow faster 
than the typical logarithmic scaling law observed in (1 + l)-dimensional quantum spin chains at 
conformally-invariant critical points [36-38]. Entanglement has thus been drastically reduced by 
the symmetry under permutations of the model, as we hinted at the beginning of the Chapter 3 . 



a One should take care with this statement, since there are other models which are symmetric under permutations 
of particles and such that the entropy of entanglement is very large, as are for instance those systems described by the 
Laughlin wavefunction [148]. 
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Figure 3.2: Entanglement entropy for N = 500 and L = 125 as a function of h and y. 



3.2.1 The y-h plane 

In order to study the different entanglement regimes, we compute the entropy in the plane spanned 
by y and h. The numerical computation can be done by taking advantage of the Hamiltonian sym- 
metries to reduce the complexity of the task to a polynomial growth in N. Results are displayed in 
Fig l3.2l for N = 500 and L = 125. For y + 1, one clearly observes a peak at the critical point h = 1 
whereas the entropy goes to zero at large h since the ground state is then a fully polarized state in 
the field direction. In the zero field limit, the entropy saturates when the size of the system increases 
and goes to S l,n = 1 for y = where the ground state approaches a GHZ-like "cat" state as in the 
Ising quantum spin chain [37,38,95, 112]. By contrast, for y = 1, the entropy increases with the size 
of the system in the region < h < 1 and jumps directly to zero at h = 1 as we shall now discuss. 



3.2.2 Analytical study of the isotropic case 

In the isotropic case (y = 1), it is possible to analytically compute the entropy of entanglement 
since, at this point, the Hamiltonian is diagonal in the Dicke basis. The ground-state energy is given 
by E (h, y = 1) = -f + jjM 2 - 2hM, with 

| WN/2), if 0<h<l 

\ N/2, if h > 1 ' y ' 



and the corresponding eigenvector is simply \N/2, M). Here, I{x) denotes the round value of x. 

To calculate the entropy, it is convenient to introduce the number n of spins "up" so that M = 
n - N/2, and to write this state in a bipartite form. Indeed, since Dicke states are completely 



3.2. Entanglement within different regimes 



41 



symmetric under any permutation of sites, it is straightforward to see that the ground state can be 
written as a sum of byproducts of Dicke states 

L 

\N/2,n-N/2) = 2^; /2 |L/2,/-L/2>® (3.5) 

1=0 

\{N - L)/2, n-l-(N- L)/2) , 

where the partition is made between two blocks of size L and (N - L) and 

L\( N-L 
I ][ n-l 

Pi = L , (3-6) 

N 

n 

defining an hypergeometric probability distribution. The expression given in Eq l3.6l corresponds 
to the Schmidt decomposition of the ground state of the system. The entropy of this state for this 
bipartition is then simply given by S^{h,y = 1) = - 2Zf=o Pi 1°§2 Pi- m tne limit N, L » 1, the 
hypergeometric distribution of the pi can be recast into a Gaussian distribution 



1 



pi^p\ = \ ^ I , (3.7) 



of mean value I ' = njj and variance 



, (N - L)L 

a- 2 = n(N - rc) ^ , (3-8) 

where we have retained the sub-leading term in (N - L) to explicitly preserve the symmetry S l,n - 
S n-l,n- The entropy then reads 

- J Alp] log 2 pf = X - (log 2 e + log 2 In + log 2 o 2 ) , (3.9) 

and only depends on its variance as expected for a Gaussian distribution b . Of course, for h > 1, 
the entanglement entropy is exactly zero since the ground state is, in this case, fully polarized in the 
magnetic field direction (n = N). For h € [0, 1) and in the limit N, L » 1, EqEH EqEHJand Eq|3~9l 
lead to 

S L , N {h,y= 1) ~ ^log 2 | L(Ar ~ L) l . (3.10) 



2 oz \ N 

Moreover, the dependence of the entropy on the magnetic field is given by 

S L , N (h, 7 = 1) - S LM (h = 0, y = 1) ~ i log 2 (l - h 2 ) , (3.11) 

and thus diverges, at fixed L and N, in the limit h — > 1~. 

b This result has also been obtained in the context of the ferromagnetic Heisenberg chain [149]. 
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N=500. L=25 
N=500. L=75 
N=500, L=125 
N=500, L=250 
N=1000. L=50 
N=1000, L=150 
N=1000, L=250 
N=1000, L=500 




Figure 3.3: Entanglement entropy at y - as a function of h for different values of N and L. Outside 
of the critical region, the entropy only depends on the ratio L/N. 



3.2.3 Numerical study of the anisotropic case 

Let us now discuss the more general situation y t 1 for which no simple analytical solution exists. 
In this case, the ground state is a superposition of Dicke states with coefficients that can be easily 
determined by exact numerical diagonalizations. Upon tracing out (N - L) spins, each Dicke state 
decomposes as in Eq l3.5l It is then easy to build the (L + 1) X (L + 1) ground state reduced density 
matrix and to compute its associated entropy. 

We have displayed in Fis l3.31 the behavior of the entropy as a function of h, for different values 
of the ratio L/N and for y = 0. For h + 1, the entropy only depends on the ratio L/N. For any y, at 
fixed L/N and in the limit h — > oo, the entropy goes to zero since the ground state becomes then fully 
polarized in the field direction. Notice that the entropy also vanishes, at h > 1, in the limit L/N — > 
where the entanglement properties become trivial. In the zero field limit, the entropy goes to a 
constant which depends on y and equals 1 at y = since the ground state is then a GHZ-like state 
made up of spins pointing in ±x directions. Close to criticality, the entropy displays a logarithmic 
divergence, which we numerically find to obey the law 



S L , N {h,y)~-a\og 2 \\-h\, (3.12) 

where a is close to 1/6 for N, L » 1 as can be seen in Fis l3.4l 

At the critical point, the entropy has a nontrivial behavior that we have studied focusing on the 
point y - which is representative of the class y ± 1. There, the entropy also scales logarithmically 
with L as in the isotropic case, but with a different prefactor. More precisely, we find 

S L , N {h = l,y * 1) ~ b\og 2 r {N ~ L) \ . (3.13) 
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N=2000, L=125 
N=2000, L=250 
N=2000, L=500 
N=2000, L=1000 



I 1 1 1 1 1 1 

1 1.02 1.04 1.06 1.08 1.1 

h 

Figure 3.4: Entanglement entropy as a function of h near the critical point for y - 0. The full line 
corresponds to the fitting law from Eg 13. 121 with a = 1/6. 



For the finite-size systems investigated here, the prefactor varies when either the ratio L/N or y is 
changed, as can be seen in Fig J3.5l However, in the thermodynamic limit Af,L » 1 (and finite L/N), 
b = 1/3 fits well our numerical results. 

In addition, at fixed L and N, the entropy also depends on the anisotropy parameter logarithmic- 
ally as 

S hN {h=\,y)-S UN {h=\,y = 0)~ f\og 2 {\-y), (3.14) 

for all -1 < y < 1 as can be seen in Fig l3.6l Here again, it is likely that, in the thermodynamic 
limit, / has a simple (rational) value which, from our data, seems to be 1/6. It is important to keep 
in mind that the limit y — > 1 and the thermodynamic limit do not commute so that Eg 13. 141 is only 
valid for y ^ 1. 

Actually, the logarithmic behavior of the laws given in Eq J3. 121 Eq 13. 131 and Eq l3.14l has been 
very recently confirmed by yet unpublished analytical computations [150], but with values of a 
and b that differ from those obtained in simulations. More precisely, it has been proven that the 
exact coefficients a and b governing the logarithmic behaviors of Eq 13 .121 and Eq.?? are 1/4 and 
1/2 respectively, instead of the values 1/6 and 1/3 obtained from the numerical computations. The 
same analytical study confirms the value of 1 /6 for coefficient / in Eq J3. 141 



3.3 Comparison to quantum spin chains 

Let us now compare the previous results with those found in the (1 + l)-dimensional quantum XY 
model. As for the Lipkin-Meshkov-Glick model, the XY quantum spin chain has two different 
universality classes depending on the anisotropy parameter. At the critical point, the entropy has 
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Figure 3.5: Entanglement entropy as a function of L at the critical point for different y and ,/V = 2000. 
The full line corresponds to the fitting law from Eq l3.13l with b = 1/3. 
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Figure 3.6: Entanglement entropy at the critical point h - 1 as a function of y. The full line 
corresponds to the fitting law from Eq l3. 141 with / = 1/6. 
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XY quantum spin chain 


Lipkin-Meshkov-Glick model 


S L (A,y = 0)~ ±log 2 (L) 
S L (A,y = 0) - S £ tf = 0,y = 0) ~ i log 2 (l - 
Si(A=l,y=l)~ ilog 2 (L) 
SlU,7= l)~4log 2 (m) 
S L {A=\,y)-S L {A=\,y=\)~\\og 2 {y) 


S L (A,y= D~ £log 2 (L) 
S L (A, r =l)-Si(/j = 0,y=l)~ilog 2 (l-A 2 ) 
Si(A=l,y = 0)~ |log 2 (L) 

5 i (ft,r = o)~-iio g2 |i-A| 
Si(* = l,y) - S £ (* = l,y = 0) ~ \ lo g2 (l - y) 



Table 3.1: Comparison of results between the XY quantum spin chain and the Lipkin-Meshkov- 
Glick model, when N » L » 1. 



been found to behave as [37, 38, 1 15] 

c ( UN-LA 
S L ,N~2 l °g2[ jf 1. (3-15) 

where c is the central charge of the corresponding (1 + l)-dimensional conformal field theory [36] 
(see Appendix B). For the isotropic case, the critical model is indeed described by a free boson 
theory with c = 1 whereas the anisotropic case corresponds to a free fermion theory with c - 1/2. 
It is striking to see that the entropy in the Lipkin-Meshkov-Glick model has the same logarithmic 
dependence with some prefactor which, as in the (1 + l)-dimensional case, only seems to depend on 
the universality class - see Eq 13 .101 and Eq J3. 131 — . Concerning the dependence with the magnetic 
field and with the anisotropy parameter, it is also worth noting that logarithmic behaviors of Eq l3.11l 
Eq l3. 121 and Eq l3.14l are similar to those found in the XY quantum spin chain [37, 38] except that 
the prefactors in the Lipkin-Meshkov-Glick model are different. A list of analogies between the 
results of the Lipkin-Meshkov-Glick model and the XY quantum spin chain in the limit N » L » 1 
is given in Table l3.ll Also, and just as a remark, it is possible to numerically check that the behavior 
of this model with respect to majorization (see Appendix A) for y + 1 and as h departs from its 
critical value is completely analogue to the case of the quantum XY model [95, 112], which was 
analytically studied in Chapter 1. Namely, the whole set of eigenvalues of the reduced density 
matrices of the ground state obey strict majorization relations as h grows, while for decreasing h 
one of the eigenvalues of the reduced density matrix in consideration drives the system towards a 
GHZ-like state in such a way that majorization is only strictly obeyed in the thermodynamic limit. 
This behavior implies a very strong sense of order of the correlations present in the ground state, in 
complete analogy to the behavior of the XY quantum spin chain. 
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3.4 Conclusions of Chapter 3 

In this Chapter we have studied the entanglement properties of a quantum spin model defined on a 
simplex. We have seen that: 

• Contrary to the intuitive idea that the quantum correlations present in the system increase 
together with the connectivity of the model, here the symmetries force the entropy to scale as 
if the system were defined on a chain. 

• Also, the Lipkin-Meshkov-Glick model presents striking similarities with the XY quantum 
spin chain: not only their phase diagrams are almost identical, but the scaling properties 
of the entanglement of the ground state seem to obey the same laws but with appropriate 
proportionality coefficients. 

The observed similarity in the behavior of this model to (1 + l)-dimensional quantum systems is 
indeed very pleasant, since quantum spin chains have been heavily studied and their properties are 
very well-known. Some of their properties seem to be directly translated into systems which, a 
priori, are not defined in (1 + 1) dimension, like the Lipkin-Meshkov-Glick model. Nevertheless, 
most of the situations that one finds when considering models which are not defined on a chain turn 
out to be much more intrincated, as we will see in the next two Chapters. Perhaps, a perturbative 
analysis around the Lipkin-Meshkov-Glick model - for instance removing a few number of links 
in the simplex and thus slightly breaking the symmetry present in the problem - could allow to 
analytically study non-trivial properties of quantum many-body systems of high dimensionality. 



Chapter 4 



Entanglement entropy in quantum 
algorithms 

The previous Chapters were focused on the properties of quantum many-body systems, basically 
from a condensed matter and field theoretical point of view. In particular, we saw that it is possible 
to apply tools from quantum information science - such as majorization and entanglement theory 
- to obtain a better understanding of the properties of these systems. We will now see that these 
tools can also be used to understand better problems arising in the area of quantum information and 
quantum computation. 

In this and the forthcoming Chapters our aim is to study a physical system which is very close 
to the spirit of quantum many -body physics: we wish to understand the properties and behavior 
of quantum computers and quantum algorithms. Indeed, a quantum computer is nothing but a 
physical system which is governed by the laws of quantum mechanics and on which we can perform 
physical actions - algorithms - such that the device is able deliver solutions to specific problems. 
Of course, the kind of problems that we can solve by using a quantum computer is necessarily 
limited by quantum physics itself, being this properly formalized by the area of quantum complexity 
theory [151]. Furthermore, it is plausible to think of a quantum computer as a device made of 
qubits which interact among themselves in some way. Therefore, a quantum computer can be 
understood as an interacting quantum many-body system. The full machinery from quantum many- 
body physics can then in principle be applied to analyze the performance of quantum algorithms. 
In particular, there is a very strong connection between quantum algorithms and quantum phase 
transitions, as we shall see. 

From the point of view of quantum computation, the design of new quantum algorithms is a 
great theoretical challenge. The most relevant property in order to understand these algorithms is 
clearly the role entanglement plays in quantum computational speedup, while some other properties 
seem to play a role as well, as we shall see in Chapter 6 with majorization [152-154]. Regarding en- 
tanglement, several results have been found [49,50, 155-159] which suggest that entanglement is at 
the heart of the power of quantum computers. An important and remarkable result was obtained by 
Vidal [49], who proved that large entanglement between the qubits of a quantum register is a neces- 
sary condition for exponential speed-up in quantum computation. To be precise, a quantum register 
such that the maximum Schmidt number of any bipartition is bounded at most by a polynomial 
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in the size of the system can be simulated efficiently by classical means. The classical simulation 
scheme proposed in [49] was, indeed, a time-dependent version of the density matrix renormaliza- 
tion group algorithm, based on the efficient updates in time of the quantum register defined in terms 
of a matrix product state [45,46]. Those methods are, indeed, tools for the classical simulation of 
the dynamics of a quantum many -body system which are also useful in the simulation of a quantum 
computation, since any quantum algorithm can be understood as the time evolution of a quantum 
many-body system [71]. Here we just sketch the basic idea of Vidal's algorithm, and leave all the 
specific details of this and other classical simulation protocols for the next Chapter. 

The figure of merit x proposed in [49] is the maximum Schmidt number of any bipartitioning 
of the quantum state or, in other words, the maximum rank of the reduced density matrices for any 
possible splitting. It can be proven that x ^ 2 S(p \ where the von Neumann entropy S(p) refers 
to the reduced density matrix of any of the two partitions. From now on, in this and also in all 
the forthcoming Chapters we shall use the following computer-science notation: the number of 
qubits in the quantum register will be denoted by n, and N - 2" denotes the dimensionality of the 
computational Hilbert space, as opposed to the condensed matter notation of the previous Chapters, 
were ./V was the number of particles present in the system. Using this notation, Vidal proved that if 
X = 0(poly(n)) at every step of the computation in a quantum algorithm, then it can be efficiently 
classically simulated. Exponential speed-up over classical computation is only possible if at some 
step along the computation x ~ exp(« a ), or S(p) ~ n b , a and b being positive constants. In order 
to exponentially accelerate the performance of classical computers any quantum algorithm must 
necessarily create an exponentially large amount of x at some point. 

As we saw in the previous Chapters, a topic of intense research concerns the behavior of entan- 
glement in systems undergoing a quantum phase transition [160]. More generally, when a splitting 
of a (d + l)-dimensional spin system is made, the von Neumann entropy of the ground state for the 
reduced density matrix of one of the subsystems S(p) - -tr(plog 2 p) at the critical point should 
typically display a universal leading scaling behavior determined by the area of the region parti- 
tioning the whole system [31, 105, 131], with at most logarithmic corrections if the system is fer- 
mionic [109-111]. As hinted in the previous Chapter, this result depends on the connectivity of the 
Hamiltonian. Using a naive reasoning, we saw there that the leading universal scaling behavior for 
the entropy of an exact bipartition of the system should typically be written in terms of the number 
of particles n as 

S(p)~«^ (4.1) 

for a (d + l)-dimensional critical non-fermionic system with sufficiently local interactions, which 
reduces to a logarithmic law for d - 1 . This explicit dependence of entanglement on dimensionality 
turns out to shed new light into some well established results from quantum computation. 

A similar situation is present in quantum adiabatic algorithms, originally introduced by Farhi et 
al. in [16], where the Hamiltonian of the system depends on a control parameter s which in turn has 
a given time dependence. The Hamiltonians related to adiabatic quantum computation for solving 
some NP-complete problems (such as 3-SAT or Exact Cover) can be directly mapped to interacting 
non-local spin systems, and therefore we can extend the study of entanglement to include this kind 
of Hamiltonians. This point of view has the additional interest of being directly connected to the 
possibility of efficient classical simulations of the quantum algorithm, by means of the protocol 
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proposed in [49]. 

Here we analyze the scaling of the entropy of entanglement in several quantum algorithms. 
More concretely, we focus on Shor's quantum factoring algorithm [8] and on a quantum algorithm 
by adiabatic evolution solving the Exact Cover NP-complete problem [16, 61-68], finding for both 
of them evidence (either analytical or numerical) of a quantum exponential speedup with linear scal- 
ing of quantum correlations - as measured by the entropy -, which seems to prohibit the possibility 
of an efficient classical simulation. We furthermore make an analytical study of the adiabatic im- 
plementation of Graver's quantum search algorithm [9,69,70], in which entanglement is a bounded 
quantity between calls to the quantum oracle even at the critical point, regardless of the size of the 
system. Let us begin, then, by considering the behavior of the factoring quantum algorithm. 

4.1 Entanglement in Shor's factoring quantum algorithm 

It is believed that the reason why Shor's quantum algorithm for factorization [8] beats so clearly 
its classical rivals is rooted in the clever use it makes of quantum entanglement. Several attempts 
have been made in order to understand the behavior of the quantum correlations present along the 
computation [157-159]. In our case, we will concentrate in the study of the scaling behavior for the 
entanglement entropy of the system. We shall first remember both Shor's original [8] and phase- 
estimation [161] proposals of the factoring algorithm and afterwards we shall move to the analytical 
study of their quantum correlations. 

4.1.1 The factoring quantum algorithm 

The interested reader is addressed to [2, 8, 161, 162] for precise details. Given an odd integer N 
to factorize, we pick up a random number a e [1,AT|. We make the assumption that a and ,/V are 
co-primes - otherwise the greatest common divisor of a and N would already be a non-trivial factor 
of N -. There exists a smaller integer r e [l,N], called the order of the modular exponentiation 
a x mod N, such that a r mod N = 1. Let us assume that the a we have chosen is such that r is 
even and a r ^ mod N + -1, which happens with very high probability, bigger than or equal to 
1 /(2 log 2 AO- This is the case of interest because then the greatest common divisor of N and a r ^ 2 ± 1 
is a non-trivial factor of N. Therefore, the original factorization problem has been reduced to the 
order-finding problem of the modular exponentiation function a x mod N, and it is at this point where 
quantum mechanics comes at work. The procedure can be casted in two different (but equivalent) 
ways: 

Shor's proposal for order-finding 

We make use of two quantum registers: a source register of k qubits such that 2 k e [N 2 , 2N 2 ], and 
a target register of n = [log 2 N~\ qubits. The quantum circuit of the quantum algorithm is shown 
in Fig l4.ll where we are making use of the Hadamard gate initially acting over the k qubits of the 
source, the unitary implementation of the modular exponentiation function 



U f \q)\x) = \q)\(x + a q ) mod AO , 



(4.2) 
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where \q) and \x) respectively belong to the source and target registers, and the quantum Fourier 
transform operator 



\m) 



(4.3) 



m=0 



All these operations can be efficiently implemented by means of one and two-qubit gates. Finally, a 
suitable classical treatment of the final measurement of this quantum algorithm provides us with r 
in few steps, and therefore the prime factorization of in a time 0((log 2 AO 3 ). 



(k) 



U 



H 



in) 



U f 



QFT 



Figure 4. 1 : Quantum circuit for the order-finding algorithm for the modular exponentiation function. 
The source and target registers have k and n qubits respectively. 



Phase-estimation proposal for order-finding 

We shall address the specific details of the generic quantum phase-estimation algorithm in Chapter 
6 and refer the interested reader to [161] for more information. For order-finding purposes, the 
quantum circuit is similar to the one shown in the previous section but slightly modified, as is 
shown in Fig |4.2l The unitary operator V/ to which the phase-estimation procedure is applied is 
defined as 

V f \x) = \(a x) mod AO (4.4) 
(notice the difference between Eq 14 .41 and Eq |4.2l ). being diagonalized by eigenvectors 

1 r_1 

|v,> = — Yj e~ 2 ™ p/r \a p mod AO (4.5) 
r P =o 

such that 

V>|v 4 > - e 2nislr \v s ) , (4.6) 

and satisfying the relation Z 4 =o ~ l^- ^ ne °P erat or is applied over the target register being 
controlled on the qubits of the source in such a way that 

A(V/)|j>|jc> - \j)V J f \x) , (4.7) 

where by A(V/) we understand the full controlled operation acting over the whole system, which 
can be efficiently implemented in terms of one and two-qubit gates. As in the previous case, the 
information provided by a final measurement of the quantum computer enables us to get the factors 
of N in a time 0((log 2 AO 3 ). 
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Figure 4.2: Phase-estimation version of the quantum circuit for the order-finding algorithm. The 
controlled operation is A(V/). The source and target registers have k and n qubits respectively. 



4.1.2 Analytical results 

We choose to study the amount of entanglement between the source and the target register in the 
two proposed quantum circuits, right after the modular exponentiation operation U / from Fig |4.1l or 
the controlled Vf operation from Fig l4.2l and before the quantum Fourier transform in both cases. 
At this step of the computation, the pure quantum state of the quantum computer is easily seen to 
be exactly the same for both quantum circuits, and is given by 

2*-l 

^ > = ^2> >|aI?modA0 ' (4 - 8) 

q=0 

and therefore the density matrix of the whole system is 

j 2*-l 

IWI = 2* Z (l?><?/|) ( la * m ° d N){aq ' m ° d Nl ) ■ (4 " 9) 

q,q'=0 

Tracing out the quantum bits corresponding to the source, we get the density matrix of the target 
register, which reads 

1 2M 

Ptarget - tr source (|^)<^|) = - J] ((p\gW\p)) modN)(a"' modN\) , (4.10) 

p,q,q'=0 

that is, 

1 2M 1 r ~ l 

Ptarget = T£ J] \a p mod N){a p mod J] \a p modN)(a p modN\ . (4.1 1) 

P =o r p=0 

The last step comes from the fact that a r mod N = I, where r e [l,N] denotes the order of the 
modular exponentiation. If 2 k were a multiple of r there would not be any approximation and the 
last equation would be exact. This is not necessarily the case, but the corrections to this expression 
are 0(1/2*), thus being exponentially small in the size of the system. 

It follows from Eq l4.1 ll that the rank of the reduced density matrix of the target register at this 
point of the computation is 

rank(ptar g et) ~ r . (4.12) 

Because r e [1,AT|, this rank is usually 0(N). If this were not the case, for example if r were 
0(log 2 AO, then the order-finding problem could be efficiently solved by a classical naive algorithm 
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and it would not be considered as classically hard. Because N is exponentially big in the number 
of qubits, we have found a particular bipartition of the system (namely, the bipartition between the 
source register and the target register) and a step in the quantum algorithm in which the entangle- 
ment, as measured by the rank of the reduced density matrix of one of the subsystems, is expo- 
nentially big. This implies in turn that Shor's quantum factoring algorithm can not be efficiently 
classically simulated by any protocol in [49] owing to the fact that at this step x — 0(N), therefore 
constituting an inherent exponential quantum speed-up based on an exponentially big amount of 
entanglement. It is worth noticing that the puipose of the entanglement between the two registers 
consists on leaving the source in the right periodic state to be processed by the quantum Four- 
ier transform. Measuring the register right after the entangling gate disentangles the two registers 
while leaving the source in a periodic state, and this effect can only be accomplished by previously 
entangling source and target. These conclusions apply both to Shor's original proposal (circuit of 
Fis l4.lt and to the phase-estimation version (circuit of Fig l4.2b . 

The behavior of the rank of the system involves that the entropy of entanglement of the reduced 
density matrix at this point will essentially scale linearly with the number of qubits, 5(p targ et) = 
log 2 r ~ log 2 N ~ n, which is the hardest of all the possible scaling laws. We will find again this 
strong behavior for the entropy in the following section, when considering an adiabatic quantum 
algorithm solving an optimization NP-complete problem. 

4.2 Entanglement in an adiabatic NP-complete optimization algorithm 

We now turn to analyze how entanglement scales for a quantum algorithm based on adiabatic evol- 
ution [16], designed to solve the Exact Cover NP-complete problem [63]. Basic background on 
NP-completeness and classical complexity theory can be found in Appendix C. We first briefly re- 
view the proposal and, then, we consider the study of the properties of the system, in particular the 
behavior of the entanglement entropy for a given bipartition of the ground state. 

4.2.1 The adiabatic quantum algorithm 

The adiabatic model of quantum optimization algorithm deals with the problem of finding the 
ground state of a given system represented by its Hamiltonian. Many relevant computational prob- 
lems, such as 3-SAT [72], can be mapped to this situation. The method is briefly summarized as 
follows: we start from a time dependent Hamiltonian of the form 

H(s{t)) - (1 - s(t))H + s{t)H P , (4. 13) 

where Ho and Hp are the initial and problem Hamiltonian respectively, and s(t) is a time-dependent 
function satisfying the boundary conditions s(0) = and s(T) = 1 for a given T. The desired 
solution to a certain problem is encoded in the ground state of Hp. The gap between the ground and 
the first excited state of the instantaneous Hamiltonian at time t will be called g(t). Let us define 
g m i n as the global minimum of g(t) for t in the interval [0, T]. If at time T the ground state is given 
by the state \Eo; T), the adiabatic theorem states that if we prepare the system in its ground state at 
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t = 0, which is assumed to be easy to prepare, and let it evolve under this Hamiltonian, then 

K£ ;7W)>| 2 > l-e 2 (4.14) 

provided that 



max 



dt 



< e (4.15) 

where H\$ is the Hamiltonian matrix element between the ground and first excited state, e « 1, 
and the maximization is taken over the whole time interval [0, T]. Because the problem Hamilto- 
nian encodes the solution of the problem in its ground state, we get the desired solution with high 
probability after a time T. A closer look at the adiabatic theorem tells us that T dramatically de- 
pends on the scaling of the inverse of g 2 min with the size of the system. More concretely, if the gap 
is only polynomially small in the number of qubits (that is to say, it scales as 0(l/poly(«)), the 
computational time is 0(poly(«)), whereas if the gap is exponentially small (0(2~ n )) the algorithm 
makes use of an exponentially big time to reach the solution. 

The explicit functional dependence of the parameter s(t) on time can be very diverse. The point 
of view we adopt in this Chapter is such that this time dependence is not taken into account, as we 
study the properties of the system as a function of s, which will be understood as the Hamiltonian 
parameter. We will in particular analyze the entanglement properties of the ground state of H(s), 
as adiabatic quantum computation assumes that the quantum state remains always close to the in- 
stantaneous ground state of the Hamiltonian all along the computation. Notice that we are dealing 
with a system which is suitable to undergo a quantum phase transition at some critical value of the 
Hamiltonian parameter in the thermodynamic limit, and therefore we expect to achieve the largest 
quantum correlations when evolving close to this point. The question is how these large quantum 
correlations scale with the size of the system when dealing with interesting problems. This is the 
starting point for the next two sections. 



4.2.2 Exact Cover 

The Exact Cover NP-complete problem is a particular case of the 3-SAT problem, and is defined 
as follows: given the n boolean variables {jc,-},-^...,,, x, - 0, 1 V i, where i is regarded as the bit 
index, we define a clause of Exact Cover involving the three qubits i, j and k (say, clause "C") by 
the equation x, + Xj + x k = 1. There are only three assignments of the set of variables [xj,Xj,x k } 
that satisfy this equation, namely, { 1 , 0, 0}, {0, 1 , 0} and {0, 0, 1 }. The clause can be more specifically 
expressed in terms of a boolean function in Conjunctive Normal Form (CNF) as 

<pc(*i, Xj, Xk) - (Xi V Xj V X k ) A (-iX, V -iX j V -iXyfc) A (~iXj V -IX j V x k ) 

A(-ocj V Xj V ->x k ) A (x t V -iXj V -a*) , (4. 16) 

so 4>c(xi, Xj, x k ) = 1 as long as the clause is properly satisfied. An instance of Exact Cover is a 
collection of clauses which involves different groups of three bits. The problem is to find a string of 
bits {x\,X2-.., x n ) which satisfies all the clauses. 
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This problem can be mapped into rinding the ground state of the Hamiltonian Hp of a spin- 1/2 
system in the following way: given a clause C define the Hamiltonian associated to this clause as 



H c = i(l+af)i(l+tr^(l+<rp 



+ I(l-af)i(l-o*)i(l-o2 
+ \{l - of)i(l + o*)l(l - a*) 

+ I(l +a f)I(l-a*)±(l-a*), (4.17) 

where we have defined o~ z |0} = |0), o~ z |l) = -|1). Note the analogy between Eg 14. 161 and Eg 14. 171 
The quantum states of the computational basis that are eigenstates of He with zero eigenvalue 
(ground states) are the ones that correspond to the bit string which satisfies C, whereas the rest of 
the computational states are penalized with an energy equal to one 3 . Now, we construct the problem 
Hamiltonian as the sum of all the Hamiltonians corresponding to all the clauses in our particular 
instance, that is to say, 

H P = H c> (4- 18 ) 

C e instance 

so the ground state of this Hamiltonian corresponds to the quantum state whose bit string satisfies 
the maximum number of clauses (all of them if the clauses are mutually compatible). We have 
reduced the original problem stated in terms of boolean logic to the hard task of finding the ground 
state of a two and three-body interactive spin Hamiltonian with local magnetic fields. Observe that 
the couplings depend on the particular instance we are dealing with, and that the spin system has 
not an a priori well defined dimensionality neither a well defined lattice topology, in contrast with 
some usual simple spin models. 

We now define our s-dependent Hamiltonian H{s) as a linear interpolation between an initial 
Hamiltonian Hq and Hp: 

H(s) = (1 - s)H + sH P (4.19) 

where we take the initial Hamiltonian Ho to be that resulting from the interaction with a magnetic 
field in the x direction: 

n , 

Ho = Yj T (1 ~ ^ ' (4 - 20) 



where dj is the number of clauses in which qubit i appears, and cr x \+) = |+), with |+) = -^=(10) + 11))> 
so the ground state of Hq is an equal superposition of all the possible computational states. Observe 



a In the next Chapter we shall consider a different implementation of H c - 
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that H(s) is, apart from a constant factor, a sum of terms involving local magnetic fields in the x and 
z direction, together with two and three-body interaction coupling terms in the z component. We 
can thus expect this system to undergo a quantum phase transition (in the limit of infinite n) as s is 
shifted from to 1. The numerical study of this phenomena is the aim of the next section. 

4.2.3 Numerical results up to 20 qubits 

We have randomly generated instances for Exact Cover with only one possible satisfying assignment 
and have constructed the corresponding problem Hamiltonians. Instances are produced by adding 
clauses randomly until there is exactly one satisfying assignment, starting over if we end up with no 
satisfying assignments. According to [63], these are believed to be the most difficult instances for 
the adiabatic algorithm. Our analysis proceeds as follows: 

Appearance of a quantum phase transition 

We have generated 300 Exact Cover instances - 300 random Hamiltonians with a non-degenerated 
ground state - and have calculated the ground state for 10, 12 and 14 qubits for different values of the 
parameter s in steps of 0.01. We then consider a particular bipartition of the system into two blocks 
of n/2 qubits, namely, the first n/2 qubits versus the rest, and have calculated the entanglement 
entropy between the two blocks. For each of the randomly generated Hamiltonians we observe 
a peak in the entanglement entropy around a critical value of the parameter s c ~ 0.7. We have 
averaged the obtained curves over the 300 instances and have obtained the plot from Fig l4.3l 

The point at which the entropy of entanglement reaches its maximum value is identified as the 
one corresponding to the critical point of a quantum phase transition in the system (in the limit of 
infinite size). This interpretation is reinforced by the observation of the typical energy eigenvalues 
of the system. For a typical instance of 10 qubits we observe that the energy gap between the ground 
state and the first excited state reaches a minimum precisely for a value of the parameter s c ~ 0.7 
(see Figi4"4l. 

We observe from Fig |4.3l that the peak in the entropy is highly asymmetric with respect to the 
parameter s. A study of the way this peak seems to diverge near the critical region seems to indicate 
that the growth of entanglement is slower at the beginning of the evolution and fits remarkably 
well a curve of the type S ~ log 2 |log 2 (s - s c )\, whereas the falling down of the peak is better 
parameterized by a power law S ~ \s—s c \~ a with a ~ 2.3, a being a certain critical exponent. These 
laws governing the critical region fit better and better the data as the number of qubits is increased. 

Analysis of different bipartitions of the system 

An explicit numerical analysis for 10 qubits tells us that all possible bipartitions for each one of the 
instances produce entropies at the critical point of the same order of magnitude - as expected from 
the non-locality of the interactions -. This is represented in Fig |4.51 where we plot the minimum 
and maximum entanglement obtained from all the possible bipartitions of the system for each one of 
the generated instances (points are sorted such that the minimum entropy monotonically increases). 

Similar conclusions follow from the data plotted in Fig l4.6l where we have considered again the 
same quantities but looking at 64 randomly-chosen bipartitions of the ground state for 10 different 
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Figure 4.3: Evolution of the entanglement entropy between the two blocks of size n/2 when a bipar- 
tition of the system is made, on average over 300 different instances with one satisfying assignment. 
A peak in the correlations appears for s c ~ 0.7 in the three cases. 
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Figure 4.4: Energies of the ground state and first excited state for a typical instance with one satis- 
fying assignment of Exact Cover in the case of 10 qubits (in dimensionless units). The energy gap 
approaches its minimum at s c ~ 0.7. 
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Figure 4.5: Minimum and maximum entropy over all possible bipartitions of a 10-qubit system for 
each of the 300 randomly generated instances of Exact Cover. Instances are sorted such that the 
minimum entanglement monotonically increases. 
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Figure 4.6: Minimum and maximum entropy over 64 bipartitions of a 16-qubit system for 10 ran- 
domly generated instances of Exact Cover. Instances are sorted such that the minimum entangle- 
ment monotonically increases. 
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Figure 4.7: Scaling of the minimum energy gap (in dimensionless units) with the size of the system, 
both in the worst case and in the mean case over all the randomly generated instances. Error bars 
give 95 per cent of confidence level for the mean. 



instances of 16 qubits. According to these results we restrict ourselves in what follows to the 
analysis of a particular bipartition of the system, namely the first n/2 qubits versus the rest. 

Scaling laws for the minimum energy gap and the entanglement entropy 

To characterize the finite-size behavior of the quantum phase transition, we have generated 300 
random instances of Exact Cover with only one satisfying assignment from 6 to 20 qubits, and 
studied the maximum von Neumann entropy for a bipartition of the system as well as the minimum 
gap, both in the worst case and in the mean case over all the randomly generated instances. We must 
point out that the scaling laws found in this section are limited to the small systems we can handle 
with our computers in an exact way. Increasing the number of qubits may lead to corrections in 
the numerical results, which should be of particular importance for a more precise time-complexity 
analysis of the adiabatic algorithm. Fig l4.7l represents the behavior of the gap in the worst and mean 
cases. From Fis l4.8l we observe that the gap seems to obey a scaling law of the type 0(l/n), n 
denoting the number of qubits, which would guarantee a polynomial-time quantum computation. 
This law is in agreement with the results in [63], and are in concordance with the idea that the 
energy gap typically vanishes as the inverse of the volume in condensed matter systems (here the 
volume is the number of qubits). Error bars in the two plots give 95 per cent of confidence level in 
the numerically calculated mean. 

We have also considered the scaling behavior of the entanglement entropy for an equally sized 
bipartition of the system, again both in the worst and in the mean case. The obtained data from 
our simulations are plotted in Fig l4.9l - where error bars give 95 per cent of confidence level in 
the mean - and seem to be in agreement with a linear scaling of entanglement as a function of the 
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Figure 4.8: Minimum energy gap (in dimensionless units) versus the inverse size of the system, both 
in the worst case and in the mean case over all the randomly generated instances. Error bars give 95 
per cent of confidence level for the mean. The behavior of the mean is apparently linear. 
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Figure 4.9: Scaling of the entanglement entropy for an equally sized bipartition of the system, both 
in the worst case and in the mean case over all the randomly generated instances. Error bars give 95 
per cent of confidence level for the mean. The data are consistent with a linear scaling. 
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size of the number of qubits. More concretely, a numerical linear fit for the mean entanglement 
entropy gives us the law S ~ O.ln. Observe that the entropy of entanglement does not saturate 
at its maximum allowed value (which would be S max = n/2 for n qubits), so we can say that 
only twenty percent of all the possible potential available entanglement appears in the quantum 
algorithm. Linearity in the scaling law would imply that this quantum computation by adiabatic 
evolution, after a suitable discretization of the continuous time dependence, could not be classically 
simulated by the protocol of [49]. Given that the scaling of the gap seems to indicate that the 
quantum computation runs in a polynomial time in the size of the system, our conclusion is that 
apparently we are in front of an exponentially fast quantum computation that seems extremely 
difficult (if not impossible) to be efficiently simulated by classical means. This could be an inherent 
quantum mechanical exponential speedup that can be understood in terms of the linear scaling of 
the entropy of entanglement. Note also the parallelism with the behavior of the entanglement found 
in Shor's algorithm in the previous section. As a remark, our numerical analysis shows that the 
quantum algorithm is difficult to simulate classically in an efficient way, which does not necessarily 
imply that the quantum computer runs exponentially faster than the classical one, as our time- 
complexity analysis is limited to 20 qubits. 

The linear behavior for the entropy with respect to the size of the system could in principle be 
expected according to the following qualitative reasoning. Naively, the entropy was expected to 
scale roughly as the area of the boundary of the splitting. This area-law is in some sense natural: 
since the entropy value is the same for both density matrices arising from the two subsystems, 
it can only be a function of their shared properties, and these are geometrically encoded in the 
area of the common boundary. For a system of n qubits, we observe again that this implies a 
scaling law for the entropy of an exact bipartition like S ~ n~ (which reduces to a logarithm for 
d - 1). Our system does not have a well defined dimensionality, but owing to the fact that there are 
many random two and three-body interactions, the effective dimensionality of the system should be 
very large. Therefore, we expect a linear (or almost linear) scaling, which is what we numerically 
obtained. While this reasoning is not valid for critical fermionic systems, it differs only by at most a 
logarithmic multiplicative correction which we did not see in our computations. The data seems to 
indicate that such an effective dimensionality is around d ~ n, thus diverging as n goes to infinity. 

It is possible to compare our seemingly linear scaling of the mean entropy of entanglement with 
the known results obtained by averaging this quantity over the entire manifold of rc-qubit pure states, 
with respect to the natural Fubini-Study measure. According to the results conjectured in [163] and 
later proved in [164], the average entropy for an equally-sized bipartition of a random «-qubit pure 
state in the large n limit can be approximated by S ~ (n/2) - 1/(2 In 2) (in our notation), therefore 
displaying as well a linear scaling law (but different from ours). In fact, this is an indicator that most 
of the «-qubit pure states are highly entangled, and that adiabatic quantum computation naturally 
brings the system close to these highly entangled regions of the pure state manifold. 

The entanglement-gap plane 

The plots in Fig l4. 101 and Fis l4.HI show the behavior of the peak in the entanglement versus the 
gap, both again in the average and the worst case for all the generated instances. Clearly, as the gap 
becomes smaller the production of entanglement in the algorithm increases. A compression of the 
energy levels correlates with high quantum correlations in the system. 
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Figure 4.10: Mean entropy of entanglement versus mean size of the energy gap (in dimensionless 
units). Error bars give 95 per cent of confidence level for the means. Each point corresponds to a 
fixed number of qubits. 
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Figure 4.12: Mean critical point for the energy gap and for the entropy. Error bars give 95 per 
cent of confidence level for the means. Note that they tend to approach as the size of the system is 
increased. 



Convergence of the critical points 

The critical point s c seems to be bounded by the values of s associated with the minimum gap and 
the maximum entropy Actually, the value of the critical point corresponding to the minimum size 
of the energy gap is systematically slightly bigger than the value of the critical point corresponding 
to the peak in the entropy. By increasing the size of the system these two points converge towards 
the same value, which would correspond to the true critical point of a system of infinite size. This 
effect is neatly observed in Fig |4. 121 which displays the values of s associated with the mean critical 
points both for the gap and for the entropy as a function of n. 

Universality 

The above results suggest that the system comes close to a quantum phase transition. The charac- 
terization we have presented based on the study of averages over instances reconstructs its universal 
behavior. Results do not depend on particular microscopic details of the Hamiltonian, such as the 
interactions shared by the spins or the strength of local magnetic fields. Any adiabatic algorithm 
solving a k-sat problem and built in the same way we have done for Exact Cover should display 
on average exactly the same properties we have found regardless of the value of k, which follows 
from universality (k - 1 is a particular case, as its Hamiltonian is non-interacting). Linear scaling of 
entanglement should therefore be a universal law for this kind of quantum algorithms. The specific 
coefficients of the scaling law for the entropy should be a function only of the connectivity of the 
system, that is on the type of clauses defining the instances. 

We have explicitly checked this assertion by numerical simulations for clauses of Exact Cover 
but involving 4 qubits (jc; + Xj + Xk + Xi - 1), which is a particular case of 4-SAT. In Fig l4.13l we 
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Figure 4.13: Entanglement as a function of the Hamiltonian parameter for clauses of Exact Cover 
involving 3 (k - 3) and 4 (k = 4) qubits, for a 10-qubit system, averaged over all the randomly 
generated instances. 



plot the behavior of the entropy of entanglement for a 10-qubit system for these type of clauses 
and compare it to the same quantity calculated previously for the clauses involving 3 qubits (the 
common Exact Cover Hamiltonian). We observe again the appearance of a peak in the entropy, 
which means that the system is evolving close to a quantum phase transition. 

Fia l4.14l and Fis l4. 151 respectively show the scaling of the energy gap in the mean and worst 
case and the scaling of the peak in the entropy in the mean and worst case as well, up to 16 qubits. 
Error bars give again 95 per cent of confidence level for the means. The behavior is similar to the 
one already found for the instances of Exact Cover involving 3 qubits (Fie l4.8l and Fig l4.9t . which 
supports the idea of the universality of the results. The minimum energy gap seems to scale in 
this case as ~ ^ (n being the number of qubits), which would guarantee again a polynomial-time 
quantum adiabatic evolution. 



4.3 Entanglement in adiabatic quantum searching algorithms 

Graver's quantum algorithm solves the problem of finding a "needle in a haystack", which is math- 
ematically defined as finding a specific element of an unsorted database by means of calls to an 
oracular function. If the database is composed of 2" elements, n being the number of bits, then the 
best classical algorithm for solving this problem takes 0(2") time as measured in calls to the oracle, 
whereas Graver's quantum algorithm takes only 0(2 n ^ 2 ) calls to the quantum implementation of the 
oracular function [9]. Optimality of Graver's quantum algorithm has been proven as well [165]. 

Let us now consider the adiabatic implementation of Graver's quantum searching algorithm in 
terms of a Hamiltonian evolution [9, 69, 70] and study its properties as a function of the number of 
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qubits and the parameter s. For this problem, it is possible to compute all the results analytically, 
so we shall get a closed expression for the scaling of entanglement. As a side remark, it is worth 
noting that the treatment made in [49] is not valid for the oracular model of quantum computation, 
as it is assumed that all quantum gates are known in advanced. Independently of this issue, we shall 
see that the system remains weakly entangled between calls to the oracle. 



4.3.1 Adiabatic quantum search 

Graver's searching algorithm [9] can be implemented in adiabatic quantum computation by means 
of the ^-dependent Hamiltonian 



H(s) - (1 - s)(I - \iff)M) + s(I - \xo)(x \) , (4.21) 

where = ^ T^Jq 1 \ x )> n * s tne number of qubits, and |xo> is the marked state. The computation 
takes the quantum state from an equal superposition of all computational states directly to the state 
|xo), as long as the evolution remains adiabatic. The time the algorithm takes to succeed depends on 
how we choose the parameterization of s in terms of time. Our aim here is to compute the amount of 
entanglement present in the register and need not deal with the explicit dependence of the parameter 
s on time and its consequences (see [69,70] for further information about this topic). 

It is straightforward to check that the Hamiltonian from Ea l4.21l has its minimum gap between 
the ground and first excited states at s - 0.5, which goes to zero exponentially fast as the number 
of qubits in the system is increased. Therefore, this Hamiltonian apparently seems to undergo a 
quantum phase transition in the limit of infinite size at s = 0.5. Quantum correlations approach 
their maximum for this value of s. 



4.3.2 Analytical results 

It can be seen (see for example [166]) that the ground state energy of the Hamiltonian given in 
Eq l4.2 ll corresponds to the expression 



EJs) = 1|1 - v /(l - 2sf + -s(] ~s) 



(4.22) 



s denoting the Hamiltonian parameter. The corresponding normalized ground state eigenvector is 
given by 



|E_(j)> = a\xo) + b ^ |jc> , (4.23) 



X^Xo 

where we have defined the quantities 
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a - a b 



b 2 = 



a - 



2" 



1 

1 + a 2 
2 n - 1 



2" 



(4.24) 



In all the forthcoming analysis we will assume that the marked state corresponds to \xo) - |0), 
which will not alter our results. The corresponding density matrix for the ground state of the whole 
system of n qubits is then given by 



Pn = b\a 2 



2a + 1)|0><0| + b 2 \(f>)((f>\ + b\a - 1)(|0)<O| + |O)<0|) , 



(4.25) 



where we have defined \<p) as the the unnormalized sum of all the computational quantum states (in- 
cluding the marked one), \<f>) = Y^=o l x ^- Taking the partial trace over half of the qubits, regardless 
of what n/2 qubits we choose, we find the reduced density matrix 



Pn/2 - b\a 2 -2a + 1)|0')<0'| + 2 n/2 b 2 \<p')(ct>'\ + b\a - 1)(|0')<O'| + |O'><0'|) , 



(4.26) 



where we understand that |0'> is the remaining marked state for the subsystem of n/2 qubits and 



\(p') = Y} x= q 1 \x) is the remaining unnormalized equally superposition of all the possible computa- 
tional states for the subsystem. Defining the quantities 



A = 



B = 



C = 



2 j_ 2 n l 2 



a 1 - + 



a 2 + 2 n 


- 1 


a + 2 n/2 


- 1 


a 2 + 2 n 


- 1 


2 n/2 




a 2 + 2 n - 


- 1 



(4.27) 



(note that A + (2 n/2 - \)C - 1), the density operator for the reduced system of n/2 qubits can be 
expressed in matrix notation as 



Pn/2 = 



A B 
B C 

{B C 



B\ 

C 

C) 



(4.28) 



in the computational basis, where its dimensions are 2 n/2 x 2 n/2 . We clearly see that the density 
matrix has a rank equal to 2. Therefore, because rank(p) > 2 5(p) Vp (where S (p) is the von Neumann 



4.3. Entanglement in adiabatic quantum searching algorithms 



67 



entropy of the density matrix p) we conclude that S(p n /2), which corresponds to our entanglement 
measure between the two blocks of qubits, is always < 1. This holds true even for non symmetric 
bipartitions of the complete system. Regardless of the number of qubits, entanglement in Graver's 
adiabatic algorithm is always a bounded quantity for any s, in contrast with the results obtained 
in the previous sections for Shor's factoring algorithm and for the Exact Cover problem. Graver's 
adiabatic quantum algorithm essentially makes use of very little entanglement between calls to the 
quantum oracle, but even this bounded quantity of quantum correlations is enough to give a square- 
root speedup. 

We have explicitly calculated the von Neumann entropy for p n ji. Because the rank of the re- 
duced density matrix is two, there are only two non-vanishing eigenvalues that contribute in the 
calculation which are 

vLt = i ( 1 ± ^1 -4(2"/ 2 - l)(AC-fi 2 )J . (4.29) 
We analyze the limit n — > oo for s + 0.5 and s = 0.5 separately. 



Entropy at s + 0.5 

In the limit of very high n we can approximate the ground state energy given in Eq l4.22l bv 



E.(s) ~ 1 (l - Vl - 4*1 - sj) . 



(4.30) 



Therefore, the quantity 



diverges at s = 0.5, which implies that this limit can not be correct for that value of the parameter. 
The closer we are to s = 0.5, the bigger is a. In this limit we find that 



a 



+ 2 n ' 2 



B 



a 2 + 2 n 

a + 2" /2 
a 2 + 2 n 



(4.32) 



(4.33) 

2 n/2 

C ~ = - , (4.34) 
a 2 + 2" 

where all these quantities tend to zero as n — > oo. It is important to note that the convergence of 
the limit depends on the value of a or, in other words, how close to s = 0.5 we are. The closer 
we are to s = 0.5, the slower is the convergence, and therefore any quantity depending on these 



68 



Chapter 4. Entanglement entropy in quantum algorithms 



parameters (such as the entropy) will converge slower to its asymptotical value. For the eigenvalues 
of the reduced density matrix we then find that when n — > oo 

A ± -> 1(1 ± 1) , (4.35) 
so A+ ~ 1 and A- ~ 0, and therefore the asymptotical entropy is 

S (s * 0.5, n -> oo) = -A+ log 2 A+ - A- log 2 A- - . (4.36) 
The convergence of this quantity is slower as we move towards s = 0.5. 



Entropy at s = 0.5 

We begin our analysis by evaluating the quantities at s = 0.5 and then taking the limit of big size of 
the system. We have that a(s = 0.5) = ^v-i ~ 2 n/2 - From here it is easy to get the approximations 



A ~ I 
2 

1 

B ~ 2^ 

C ~ — =-r, (4-37) 

2«/2+i v ' 



and therefore 



1 

so A+ — > ^ and S(s = 0.5, n — > oo) = 1. According to Eg 14.381 we can evaluate the finite size 
corrections to this behavior and find the scaling of the entropy with the size of the system for very 
large n. The final result for the entropy at the critical point reads 

S(s = 0.5,72 » 1)~ l--Ar2-" /2 . (4.39) 

In 2 

Note that the entropy remains bounded and tends to 1 for s = 0.5 as a square root in the exponential 
of the size of the system, which is the typical factor in Graver's quantum algorithm. 

We have represented the evolution of the entanglement entropy as a function of s for different 
sizes of the system in Fig |4.16l and have plotted in Fig l4. 17l the maximum value of the entropy along 
the computation as a function of the size of the system according to the expression given in Eg 14.391 
We can now compare the two plots with Fis l4.3l and Fig l4.9l in the previous section. The behavior for 
the entropy in Graver's adiabatic algorithm is dramatically different to the one observed in the NP- 
complete problem. Entanglement gets saturated in Graver's adiabatic algorithm even at the point 



1 + Jl -4 2" J |i— - — 

\ 42 n/2 2 " 
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(4.38) 
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Figure 4.16: Von Neumann entropy for the reduced system as a function of s for 10, 12 and 14 
qubits. As the size of the system increases the entropy tends to zero at all points, except at s = 0.5 
in which tends to 1 . 
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Figure 4.17: Von Neumann entropy for the reduced system at s - 0.5 as a function of n. For infinite 
size of the system there is a saturation at 1 . 
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at which the gap vanishes, which reminds us of short ranged quantum correlations in non-critical 
quantum spin chains' 3 . 

Let us note that, in the limit of infinite size, the quantum state in Graver's algorithm is separable 
with respect to any bipartition of the system (and therefore not entangled, as it is a pure state) for 
any s except for s = 0.5. All the entanglement along the algorithm is concentrated at this point, but 
this entanglement is still a bounded quantity and actually equal to 1 . Consequently, a small amount 
of entanglement appears essentially only at one point when the size of the system is big, whereas 
the rest of the algorithm needs to handle just separable states. We point out that these results apply 
as well to the traditional discrete-time implementation of Graver's searching algorithm, as the states 
between iterations are the same as in the adiabatic version for discrete s values. 

4.4 Conclusions of Chapter 4 

In this Chapter we have studied the scaling of the entanglement entropy in several quantum al- 
gorithms. In order to be precise: 

• We have analytically proven that Shor's factoring quantum algorithm makes use of an expo- 
nentially large amount of entanglement in the size of the system between the target register 
and the source register after the modular exponentiation operation, which in turn implies the 
impossibility of an efficient classical simulation by means of the protocol of Vidal in [49]. 

• We have provided numerical evidence for a universal linear scaling of the entropy with the 
size of the system together with a polynomially small gap in a quantum algorithm by adia- 
batic evolution devised to solve the NP-complete Exact Cover problem, therefore obtaining 
a polynomial-time quantum algorithm which would involve exponential resources if simu- 
lated classically, in analogy to Shor's algorithm. Universality of this result follows from the 
fact that the quantum adiabatic algorithm evolves close to a quantum phase transition and 
the properties at the critical region do not depend on particular details of the microscopic 
Hamiltonian (instance) such as interactions among the spins or local magnetic fields. 

• We have also proven that the von Neumann entropy remains bounded by 1 between calls to 
the quantum oracle in Graver's adiabatic algorithm regardless of the size of the system and 
even at the critical point. More concretely, the maximum entropy approaches one as a square 
root in the size of the system, which is the typical Graver's scaling factor. 

Our results show that studying the scaling of the entropy is a useful way of analyzing entangle- 
ment production in quantum computers. Results from the study of quantum many -body systems can 
be directly applied to bring further insight into the analysis of the quantum correlations present in a 
quantum computer. Different entanglement scaling laws follow from different situations according 
to the amount of correlations involved, as can be seen in Table 14.1 1 A quantum algorithm can be un- 
derstood as the simulation of a system evolving close to a quantum phase transition. The amount of 

b A somehow similar situation is present in (1 + l)-dimensional quantum spin chains outside of the critical region, 
where the entanglement entropy also reaches a saturation when increasing the size of the system [22,37,38]. Saturation 
does not appear in higher dimensional systems. 



4.4. Conclusions of Chapter 4 



71 



entanglement involved depends on the effective dimensionality of the system, which in turn governs 
the possibilities of certain efficient classical simulation protocols. 



Problem 


Scaling of the entanglement entropy 


Adiabatic Exact Cover's quantum algorithm 


S = 0(n) 


Shor's quantum factoring algorithm 


S = 0(log 2 r) ~ 0(n) 


Critical (d + l)-dimensional fermionic lattices 


S - 0{wN~ log 2 n) 


Critical (d + l)-dimensional bosonic lattices 


S = Oin*?) 


Critical (1 + l)-dimensional spin chains 


S = 0(log 2 n) 


Non-critical (1 + l)-dimensional spin chains 


S = 0(1) 


Adiabatic Graver's quantum algorithm 


S = 0(1) 



a 

ID 

£ 

"bo 



Table 4.1: Entanglement scaling laws in different problems, in decreasing complexity order. 



These scaling laws provide also a new way of understanding some aspects from one-way quantum 
computation. It is known that the so-called cluster state of the one-way quantum computer can be 
generated by using Ising-like interactions on a planar (2 + l)-dimensional lattice [167-169]. This 
fact can be related to the at least linear (in the size of a box) behavior of the entropy for spin sys- 
tems in (2 + 1) dimensions. (1 + l)-dimensional models seem not to be able to efficiently create 
the highly-entangled cluster state. Again, this fact can be traced to the logarithmic scaling law of 
the entropy in spin chains which is insufficient to handle the large amount of entanglement to carry 
out for instance Shor's algorithm. Note also that (d + l)-dimensional systems with d > 3 bring 
unnecessarily large entanglement. 

Quantum phase transitions stand as demanding systems in terms of entanglement. They are very 
hard to simulate classically. It is then reasonable to try to bring NP-complete problems to a quantum 
phase transition setup, which quantum mechanics handles naturally. 



Chapter 5 



Classical simulation of quantum 
algorithms using matrix product states 

In Chapter 4 we saw that understanding the detailed behavior and properties of quantum many- 
body systems plays a role in different areas of physics. Those systems whose properties can be 
analytically found are typically called integrable systems and offer a way to study, for instance, the 
low-energy sector of different models. It is a pity, though, that many of the models that we know 
are not integrable, in the sense that it is not even known whether it is possible or not to study in an 
exact way their fundamental properties. The realistic alternative is, then, to use different techniques 
based on numerical simulations by means of computer programs, so that we can get a detailed 
understanding of the system. 

While it is possible to numerically study the low-energy properties of any model by means of 
an exact diagonalization of the quantum Hamiltonian or related techniques, such a possibility is 
always limited to a relatively small number of particles due to the exponential growth in the size 
of the Hilbert space. Indeed, this is at the heart of the motivation to build a quantum computer, as 
originally proposed by Feynman [1]. Using standard present technology, a faithful numerical study 
of the ground-state properties of a general quantum Hamiltonian can be achieved for systems up to 
the order of 20 spins, as we did in the previous Chapter. Luckily enough, other numerical techniques 
are possible. For instance, quantum Montecarlo algorithms have provided good results for some 
systems while they fail for some others due to the presence of the so-called sign problem [170]. 
Another example of successful numerical technique has been the density matrix renormalization 
group (DMRG) algorithm, as introduced by White in [20]. While it was soon realized that DMRG 
produced extremely accurate results when computing the ground-state energy of quantum systems 
in one spatial dimension, it was also realized that the method did not work so well when applied to 
higher dimensional systems [171, 172]. Even in the (1 + l)-dimensional case, there was a difference 
in the performance of the algorithm between open and periodic boundary conditions, and between 
non-critical and critical systems, the former being the more successful in both cases. Nevertheless, 
DMRG has been the algorithm of reference for computing the low-energy properties of quantum 
models with one spatial dimension during the last decade. 

After the appearance of DMRG, a notorious result was found by Ostlund and Rommer in [47], 
where they showed that the original DMRG algorithm can be completely understood in terms of 
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the so-called matrix product states. Originally introduced in the valence-bond model of Affleck, 
Kennedy, Lieb and Tasaki [45,46], generalized by Fannes, Nachtergaele and Werner [48], and re- 
discovered in the field of quantum information science by Vidal [49], matrix product states have 
been proved to be an extremely useful tool in order to develop numerical techniques for comput- 
ing the low-energy properties together with the dynamics of sufficiently local Hamiltonians in one 
spatial dimension [50-57], and have inspired as well several numerical techniques to study higher- 
dimensional systems [58-60]. 

The natural question arises of whether matrix product states can be applied to simulate the 
dynamics of a quantum computer. The content of this Chapter is aimed to show that this is indeed 
possible and that we can handle relatively large simulations with controlled accuracy. We call the 
parameter controlling the size of the matrices x, which was already introduced in Chapter 4, and 
which can in turn be related to the entanglement entropy S of a considered bipartition of the system 
like* >2 S{ p\ As we shall see, the total time cost of the simulation scales polynomially in parameter 
X- Thus, we expect this approximation scheme to fail whenever the inherently needed x is 0(2 n ), 
n being the number of qubits of the quantum register. Nevertheless, it may be possible in some 
of these cases that by keeping only x - O(po\y(n)) in the simulation we already get a reasonable 
approximation to the exact computation. This is indeed the case of the quantum algorithm that we 
consider here. We study the numerical performance of the classical simulation scheme for quantum 
computations originally proposed by Vidal in [49] based on matrix product states, when applied to 
the simulation of an adiabatic quantum algorithm solving the Exact Cover NP-complete problem. 
The performance of this quantum algorithm was already addressed in Chapter 4, where we saw 
that the typical entanglement entropy of the system for a given bipartition tends to scale roughly 
as S ~ O.ln, which makes the parameter* exponentially big in the number of qubits and thus 
forbids the possibility of an exact classical simulation. Nevertheless, the fact that the coefficient in 
front of the scaling law of the entropy is small inspires us to think that, perhaps, it should still be 
possible to perform a relatively good approximated classical simulation of the quantum algorithm 
by keeping a small amount of x along the evolution. Notice that this is a necessary, while not 
sufficient condition to have a good approximation of the evolution of the quantum algorithm. Let 
us then proceed in what follows with an explanation of what matrix product states are and how do 
they inspire numerical simulation algorithms for time evolution, moving then to our explicit results 
on the numerical simulation of a quantum computer. 



5.1 The matrix product state ansatz 

A matrix product state is a parameterization of a pure quantum state of n local systems (like, for 
instance, qubits) in terms of the amount of bipartite entanglement present in the state. Here we 
derive this ansatz from two different perspectives: on the one hand, we show how matrix product 
states appear from the point of view inspired by Affleck, Kennedy, Lieb and Tasaki in [45,46] based 
on projectors on some ancillary unphysical particles; on the other hand, we show how it is possible 
to obtain a matrix product state by means of a series of Schmidt decompositions of the quantum 
state at hand, in the way done by Vidal in [49]. These two perspectives complement each other, and 
give different insights about the significance of the different parameters and quantities that appear 
in the ansatz. 
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Derivation by means of projectors 

Let us consider a set of n physical local J-level systems, described by (pure) quantum state given by 

d d d 

W) = ^ 2u ' ' ' Zj c hh,~,i n \h,h, •■-,'«>, (5.1) 

i! =1(2 = 1 J„=l 

where the states |//), / = 1,2, ... ,n denote a local J-level basis, and ci u i 2 ^j n are the corresponding 
d n coefficients specifying the state. We now consider the following picture. First, imagine that the 
local systems are placed on a linear chain. Second, let us represent the physical local <i-level systems 
by means of two ancillary unphysical particles, each one of them being described by a Hilbert space 
of dimension^, together with a projector from the joint ancillary Hilbert space of dimension^ 2 to 
the physical Hilbert space of dimension d. We also assume that the state of the ancillary particles 
(without the projectors) is in a dimerized state of maximally entangled pairs of dimension x- The 
projector on the local Hilbert space at site / is represented by the three-index tensor 

A Wi (C ^ 

where the index ij = 1, 2, . . . , d corresponds to the physical local Hilbert space, while the indexes 
ai-i = 1,2, ,..,x and a\ = \,2,...,x correspond to the two ancillary Hilbert spaces. This is 
represented in Fig l5.ll 

projection on physical local d-level system 




Figure 5. 1: Graphical representation of a matrix product state in terms of projections. The projectors 
act on a dimerized state of maximally entangled pairs of dimension x- 

At every site, and for each value of the physical index, we have then a matrix. Because the 
ancillary particles are in a dimerized state of maximally entangled pairs, the coefficients c, M - 2 ,...,,„ of 
the system are decomposed as products of matrices, hence the name of matrix product state. The 
explicit form of the state is 

W {a} 

where the sums are to be understood from now on over the complete range of the set of physical 
indices [i] and ancillary indices {a}. 
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Derivation by means of Schmidt decompositions 

Consider again the same set of n physical local J-level systems described by the pure state of Eq l5.ll 
where we assume as before that the local systems are sorted from 1 to n in such a way that they can be 
thought as placed on a linear chain. Following Vidal [49], if we perform the Schmidt decomposition 
between the local system 1 and the remaining n - 1 we can write the state as 

lf> = £4W>l<"' n) >, (5-4) 

where are the Schmidt coefficients, \<p„i) an d \<p®i w& me corresponding left and right 
Schmidt vectors, and ct\ = 1, 2, . , . ,d. By expressing the left Schmidt vector in terms of the original 
local basis for system 1 the state can then be written as 

l^^^^'l, (5.5) 

km 

r^J" being the appropriate coefficients of the change of basis, that is, \<p„^) = X (1 ^la^l'i)- At this 
point, we expand each Schmidt vector \(f> a [ ) in the original local basis for system 2, that is, 

l<' n) > = ^ fe>l<^> - (5.6) 

h 

We now write the unnormalized quantum state \^J"^) in terms of the at most d 2 eigenvectors of the 
joint reduced density matrix for systems (3,4,..., n), that is, in terms of the right Schmidt vectors 
\<Pa 2 ) of the particular bipartition between the first two local systems and the rest, together with 
the corresponding Schmidt coefficients A^: 

Replacing the last two expressions into Eq l5.5l we obtain 

By iterating the above procedure we finally get a representation of the quantum state in terms of 
some tensors T and some vectors A: 

I, /,\ _ V V p(l)!l 3 (l) r (2)i2 i(2) i(n-l) r (n)i n ,. , , > r - Q x 

{/} la) 

Several remarks are to be considered at this point. First, notice that the above decomposition 
immediately provides the Schmidt vectors A of all the possible contiguous bipartitions of the system. 
Second, the state from Eq |5.9l is indeed a reparametrization of a matrix product state of the form 
given in Eq l5.3l if we define the matrices at site / in the following way: 



m = r m 1(0 

(Xi-ioci ~ 1 ai-.\ai"~ai 



(5.10) 
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Third, we see that the maximum allowed rank of the different indices ai, I = 1, 2, . . . ,n — 1, is 
site-dependent, since the size of the Hilbert spaces considered when performing the consecutive 
Schmidt decompositions depends on the site. In particular, we have that, at most, ai = 1,2, ... ,d' 
for I = 1,2,..., L«/2J, and a, = 1,2,..., d (n ' l) for / = [n/2\ + 1, |n/2J + 2, ...,n- 1. Actually, 
the fact that the maximum allowed range of the matrix indices is site-dependent can also be seen 
from Eq l5.3l bv performing an appropriate set of concatenated singular value decompositions of the 
matrices defining the state. In practice, however, many of the Schmidt coefficients for the different 
contiguous bipartitions of the system shall be equal to (or almost equal to) zero depending on the 
particular state being considered. Let us then call xiU P) the local Schmidt rank for the bipartition 
between the / and the / + 1 sides for a given permutation P of the particles. We shall now define^ 
as the maximum Schmidt rank over all the possible bipartitions of the system, that is 

* = max *(/,*>)• (5.11) 

We immediately see from this definition that the parameter x controlling the maximum possible 
size of the matrices in a matrix product state of n particles is, indeed, a measure of the maximum 
bipartite entanglement that is present in the system. This representation is very appealing, since 
quantum states with low (bipartite) entanglement can then be represented by small matrices, while 
highly-entangled states must necessarily be described by matrices of large size, corresponding to 
the idea that the more entangled a system is, the harder it is to perform an exact classical description 
of it. 

The above picture can be made specific by noticing that x ^ d s , where S is the entanglement 
entropy (measured in e-dits) corresponding to any possible bipartition of the system. The study of 
the scaling of the entanglement entropy can thus be translated into the study of the possibility or not 
of an efficient representation of the quantum state in terms of a matrix product state. To be more 
precise, matrix product states allow a representation of the state in terms of 0(ndx 2 ) parameters 
instead of the original d n coefficients. Therefore, those quantum states with^f = 0(poly(«)) can be 
efficiently classically represented by a matrix product state, while those where x - 0(2") cannot. 
In fact, the computation of the expected values of local observables can be done in 0(x 3 ) time, 
thus being efficient for systems with small enough x- This is an important property, since it means 
that the matrix product state representation is not only nice, but useful as well, in the sense that it 
allows to compute important physical quantities, like correlators, in an efficient way. Any possible 
parameterization of a quantum state which does not allow to efficiently compute physical properties 
is not a useful parameterization for computational purposes. How to efficiently compute correlators 
with matrix product states can be found for instance in [57]. 

The matrix product state parameterization has been very useful in computing low-energy prop- 
erties of some sufficiently local Hamiltonians, and also the dynamics of quantum states. We shall 
not explain here the details of some optimization algorithms like DMRG, and the interested reader 
is addressed to the huge amount of existing literature about this (see for example [57, 173]). We 
do sketch, however, the basic ideas on how to proceed for computing dynamical evolutions with 
matrix product states. In fact, some optimization algorithms, like euclidean time evolution, can also 
be understood in terms of the dynamical procedures that we explain in what follows. 
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5.1.1 Computing dynamics 

In this section we explain how to compute the dynamics of a matrix product state. Our model for 
dynamical evolution is based on the application of a set of unitary gates acting either on one or two 
local J-level systems, which could perfectly correspond to a discretization of the continuous time 
evolution driven by a generic one and two-body Hamiltonian. 

Let us begin this explanation by considering the effect of a unitary gate U {t) acting over a single 
J-level system /. The consequence of this operation involves an updating of the matrices A (/) at site 
/ that goes as follows: 

A'w-iai = ^ Uf^ai-m • (5-12) 
h 

Notice that this type of local gates does not affect the ancillary indices. Entanglement is thus 
unaffected, which is a necessary condition since we are just performing a local operation. 

The effect of non-local unitary gates acting on different local systems is less obvious. We ini- 
tially consider the case of a non-local gate £/ (7 ' /+1) involving contiguous local systems / and I + 1. 
Let us define 

Z TJ Q,M) Al)h _ fJ'/i+i (C io\ 

a '- iai a ' a ' +i ~ - yJ-iJ) 

Unlike with local gates, the action of an interacting gate does not preserve the product form 
of the tensors A. To reestablish the matrix product state structure we need to rewrite using a 
Schmidt decomposition. The procedure to follow is to compute the reduced density matrix from the 
bipartition of the system between the / and I + 1 sides, which for the I + 1 side reads 

n ij = V U (I ~ l) \ 2 ® ki ®* kj (5 14) 

k,ai-i 

where we have made use of the at most x known Schmidt coefficients rf&jj? for the cut between 
the / - 1 and the / sides. After diagonalizing p using (i,ai+\) and 0',/?/+i) as composed indices, 
we directly read from the eigenvalues the at most dx updated Schmidt coefficients X® for this 
bipartition, together with the updated matrices A'^^'* 1 from the coefficients of the eigenvectors. 

Finally, the new tensors for system / are easily calculated as A'^ m - ^'a,a!+i !+1 ®a'_iWr 
Non-local gates between non-contiguous systems can be reduced to the previous case by using 
SWAP gates, producing a typical overhead of 0(n) operations. Notice that all our manipulations 
can be done in a time that grows like 0(x 3 ). 

As we have seen, non-local gates entangle the system by increasing the size of the matrices that 
must be kept in the classical simulation scheme. Each time an entangling gate is operated on two 
neighboring systems, the index of the connected ancillae is multiplied by d. To keep the numerical 
simulation under control, a (non-unique) truncation scheme is needed to stop the exponential growth 
of ancillary dimensions. The ability in this truncation is the key element for the success of the time- 
evolution algorithm. Here we explain two possible truncation schemes, the first one based on the 
original proposal of Vidal [49] of an optimal local truncation, and the second one inspired on the 
methods of Verstraete and Cirac [53-55, 58] based on an optimal non-local truncation procedure. 
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Before entering into the details of the possible truncation schemes, let us introduce a graphical 
representation of the quantum state that shall be useful in what follows. We represent the tensor 
Aa^a, at site / by the following diagram: 

a l 

h 

Figure 5.2: Diagrammatic representation of the tensor at site I. 

With this notation, a matrix product state like the one from Eq l5.3l is represented by means of 
the following tensor network: 

t f f t ? -t— f 

Figure 5.3: Diagrammatic representation of a matrix product state in terms of a tensor network. 

In the above figure we have decided to drop off the name of the indices of the matrices since 
they do not bring any extra information. Each one of the dots represents a specific particle. Vertical 
lines correspond to the indices of the physical Hilbert spaces and run up to d, while horizontal links 
between the dots correspond to the ancillary indices and run at most up to x- Now, we are in a 
position to discuss the different truncation procedures. 

Local truncation scheme 

After the application of a non-local gate on the adjacent systems I and I + 1, the obtained matrix 
product state is identical to the original one with the only exception that matrices for sites / and / + 1 
have been updated, and the rank of the link connecting these two matrices has been multiplied by 
d. A possible truncation procedure is to only change the matrices at sites / and / + 1 , computing 
two new matrices with ancillary indices up to x, in such a way that the difference with the original 
state is minimum (or analogously, the overlap with the original state is maximum). This is a local 
scheme, since it only affects the two very specific matrices of the whole matrix product state that 
were touched by the action of the unitary gate. It is easy to see that optimality in this truncation is 
achieved by keeping the x terms in the range of the common index that correspond to the largest 
eigenvalues l/l'^l 2 of the reduced density matrices of the bipartition of the system between the sites 
/ and l+l. The diagrammatic representation of this truncation is shown in Fig l5.4l 

Notice that given the locality of the procedure, this scheme seems to be a good way to implement 
a truncation in order to eventually parallelize the code of the classical simulation algorithm. More 
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/ l+l 

t f f 9-*?~ - "t— ? 



t 

/ l+l 



ft? 




Figure 5.4: Local truncation scheme. Black dots correspond to old matrices, white dots correspond 
to updated matrices after the unitary evolution, and the thick link line has a rank at most d\- Only 
matrices at sites / and / + 1 are truncated (indicated by dashed dots), and this is done by keeping 
only the most relevant x terms of the corresponding Schmidt decomposition. 



precisely, one could think of different nodes of a computer network, each one of them storing one 
matrix (or a finite set of them). This truncation scheme would only involve communication between 
the two nodes on which the non-local gate operates, leaving the rest n — 2 nodes untouched, and 
therefore involving a small amount of information to be sent between different nodes. 



Non-local truncation scheme 

Given the above local procedure, we can see that there exists the possibility to improve the precision 
in the truncation by means of a non-local updating of the matrices that define the matrix product 
state. The main idea is as follows: instead of performing an optimal truncation only in matrices at 
sites / and l+l, perform an optimal truncation in all the matrices defining the matrix product state, 
that is, find a new state with new matrices for all the sites with ancillary indices up to % such that 
the distance to the original state is minimum. This is represented in Fig l5.5l 

In order to find the new optimal matrices it is possible to proceed in the following way. Let us 
call \ifj') the exact state after the non-local unitary evolution, and |i/r) the new matrix product state 
that we use to approximate We wish to maximize the quantity (^|i/r)| over all possible matrix 
product states with ancillary indices up to^ with the normalization constraint K^l^)! 2 = 1. In 
order to perform this minimization, we fix all the matrices of |i/r) to a fixed value except the first 
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/ l+l 

f t t ?— ?~ ^' -f—f 



1 

I l+l 




Figure 5.5: Non-local truncation scheme. Black dots correspond to old matrices, white dots corres- 
pond to updated matrices after the unitary evolution, and the thick link line has a range dx- We find 
new matrices at every site (indicated by dashed dots) with ancillary indices up to x sucn that the 
distance to the original state is minimized. 



one, and maximize the overlap with respect to the first matrix with the appropriate normalization 
constraint, which can be done in 0(x 3 ) time 3 . Once the values of the first matrix are found, we re- 
peat the procedure maximizing with respect to the second matrix and finding a better approximation 
to the original exact state. The complete maximization is then performed by repeating this proced- 
ure sequentially for every site, and sweeping back and forth along the system until some desired 
convergence is achieved. 

Indeed, this truncation scheme does not require the non-local gate to be necessarily applied on 
adjacent systems. Imagine that we wish to apply a non-local gate £/ (/ ' m) between distant systems 
/ and m. It is possible to see that any such unitary matrix J/^ m) e \J(d 2 ) can always be written 
as U^ l,m) - Y^a,b CabO^a ® 0^ n \ where 0® and are 2d 2 local operators acting respectively 
on sites / and m (d 2 operators per site), and C a b are d 4 coefficients 13 . Performing a singular value 
decomposition of the coefficient C a b, this can be written as C a b - Y^i U a ^D^V M b, and therefore 
the original unitary matrix can be expressed as £/ (, ' m) - £ O® ® where we have defined 



a This is valid for the case of open boundary conditions that we analyze here. Periodic boundary conditions may 
involve a larger computational time than our case. 

b It is possible to see this property by expressing the unitary operator as the exponential of a local basis for the algebra 
u(d 2 ) and performing a Taylor expansion. 
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the operators O® = £ a Ua^O^D 1 ^ 2 and - Yib V^O^D^ 2 . Applying these operators on the 
original matrix product state is equivalent to redefine the tensors at sites / and m in such a way that 
we add a new index /i of rank d 2 : 

Z^fl?;?* (5-15) 

u 

y pim)hn Q(m) 
Z—l ff,„-ia„, fj ■ i' m i m 

•m 

Given the above equation, we see that after the application of the unitary gate, the sites / and m get 
linked by a common index fi. This is another way of understanding how non-local gates entangle 
the system, namely, by creating new bonds between the sites on which they act. Now, it is possible 
to perform again a non-local truncation much in the same way as before, by finding new matrices 
for all the sites with only two ancillary indices up to x an d also in C?Cf 3 ) time as well. This is 
represented in Fie l5.6l 




Figure 5.6: Non-local truncation scheme. Black dots correspond to old matrices, and white dots 
correspond to updated tensors after the unitary evolution. The action of a non-local gate has created 
a new link of rank d 2 between sites / and m. We find new matrices at every site (indicated by dashed 
dots) with ancillary indices up to x suc h that the distance to the original state is minimized. 

We shall expect better accuracies for this non-local truncation scheme than for the local trunca- 
tion procedure, basically because we optimize over a larger set of parameters, and because we do not 



fi ; <xi-\ai 
1 P ; a m -\a m 
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have to necessarily implement SWAP operations in order to perform non-local gates between distant 
systems, thus reducing the number of truncations to be applied in the simulation. Nevertheless, this 
scheme has the drawback that the number of operations to be done at each truncation step is bigger 
than in the local case. Also, the fact that the truncation is non-local makes it a bad candidate for a 
possible parallelization of the numerical code, since all the nodes of the computer network should 
communicate among themselves at each truncation step in order to perform the approximation of 
the exact state by a new matrix product state. 

5.2 Classical simulation of an adiabatic quantum algorithm solving 
Exact Cover 

In this section we show the results of a simulation of a quantum algorithm using matrix product 
states. More precisely, we have implemented the local truncation scheme explained in the previous 
section to the simulation of a quantum adiabatic algorithm solving hard instances of the Exact 
Cover NP-complete problem. The performance of this algorithm was already analyzed in detail 
in the previous Chapter by means of an exact numerical computation of its properties up to 20 
qubits. There we saw that the entanglement entropy of a typical bipartition of the system seems 
to scale as S ~ O.Ira, ra being the number of qubits. We also found that the linear scaling of the 
entanglement entropy forbids the possibility of an efficient numerical simulation with the methods 
of [49]. The reason becomes clear now, since a linear scaling of the entanglement entropy involves 
an exponentially big^ in the number of qubits, and therefore any algorithm based on matrix product 
states must necessarily handle matrices of exponential size in order to get a result sufficiently close 
to the exact one. In any case, the possibility of a numerical simulation of this quantum algorithm by 
using matrix product states is motivated in part by the fact that the coefficient of the scaling law for 
the entanglement entropy seems to be rather small (only 0.1). Thus, even though we should need 
an exponentially big x to perform a very accurate simulation of the adiabatic quantum algorithm, it 
could be possible that already good simulations can be performed by keeping a relatively small x- 
Furthermore, the performance of a classical simulation of a quantum algorithm by using the matrix 
product state ansatz may bring further insight on the way entanglement is used along the quantum 
evolution. As we shall see, the basic features of the quantum algorithm can still be observed even 
in the case of a highly-truncated simulation with very small x- 

Let us sketch the basic features of our simulation. First, let us remind that classically hard 
instances of Exact Cover seem to appear at the so-called easy-hard-easy transition around m ~ 
0.8ra [174], m being the number of clauses and ra being the number of qubits. We have generated 
such hard instances, with the additional property of having only a unique satisfying assignment. 
The generation of hard instances is in itself a hard problem for which we have developed specific 
algorithms, essentially based on the iterative addition of random and non-redundant clauses until 
the number of solutions of the instance is one. The quantum algorithm for a given Exact Cover 
instance follows the adiabatic evolution of the ground state of a Hamiltonian defined by H(s) = 
(1 - s)Hq + sHp, where the adiabatic parameter is s = t/T and t runs up to a total predetermined 
time T. We take the initial Hamiltonian to be Hq = £" =1 y(l -erf) where di stands for the number of 
clauses where qubit i enters. The non-local problem Hamiltonian corresponds to the sum of clauses 
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defined as 

H P = 2 (Zi + Zj + Zk-l) 2 , (5.16) 

C 6 instance 

where z, ■ = (1 - cr?)/2 has eigenvalues and 1, and C stands for a clause involving bits i, j and k 
in the specific instance. Notice the difference between the problem Hamiltonians from Eq 15 .161 and 
from Eq l4.18l Both Hamiltonians describe correctly the solution to an Exact Cover instance in its 
ground state. The essential difference between them is that while the Hamiltonian of Eq l4. 1 81 has 
three-body interactions, the Hamiltonian of Eq l5.16l has not. The problem Hamiltonian that we use 
in this Chapter is built only from one and two-body terms, together with local magnetic fields, and 
its evolution can therefore be classically simulated by the algorithms based on matrix product states 
that we have already discussed, based on the efficient updatings of the register after performing one 
and two-body unitary gates. At the level of eigenvalues, notice that the only difference between 
the two Hamiltonians comes on the eigenvalues of the excited states, thus keeping the properties of 
the low-energy sector untouched. In fact, it is easy to see by means of direct simulations that an 
adiabatic quantum algorithm based on this problem Hamiltonian shows the same important features 
as the ones already described in Chapter 4, in particular the appearance of a quantum phase transition 
at s c ~ 0.69 in the thermodynamic limit, together with a linear scaling of the entanglement entropy 
with the number of qubits with a small scaling coefficient of the order of 0.1. 

Exact simulations of quantum algorithms by adiabatic evolution solving hard instances of sat- 
isfiability problems were carried so far up to 30 qubits [175]. Here we present the possibility of 
performing approximated simulations of this quantum algorithm beyond that number. 

5.2.1 Discretization of the continuous time evolution in unitary gates 

Let us now turn to discuss the detailed way matrix product states can handle the simulation of the 
adiabatic evolution of Exact Cover. The simulation needs to follow a time evolution controlled by 
the s-dependent Hamiltonian. This continuous unitary time evolution from time to time T can be 
discretized as follows: 

U(T, 0) = U(T, T - A) . . . U(2A, A)U(A, 0) , (5.17) 

where the increment A = £ defines the discretization, M being a positive integer. Our simulations 
indicate that we can take the value A = 0.125 while keeping sufficient accuracy - as compared to 
smaller A - in all of them. After / steps s = j = y = jj, that is I = 0, . . . M. 

At any point / along the evolution the unitary operator U((l + 1)A, /A) needs further subdivision 
into elementary one and two-qubit gates. This requires the use of Trotter's formula to second order 
[176-178]: 

U((l + 1)A, ZA) = e iAH(s) ~ {jkQ-M>^*Hpjha-M>y ? (5 .1 8) 

where the partition in Hq : Hp : Hq minimizes the number of two-qubit gates as compared to the 
alternative partition Hp : Hq : Hp. We have verified as well that we can maintain a faithful classical 
simulation by choosing 6 = A. Notice that the split of exponentials in the Trotter's expansion of 
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Eq l5.18l is chosen so that Ho is explicitly separated from Hp, so that this brings the advantage that 
both pieces of the Hamiltonian can be decomposed in mutually commuting one and two-qubit gates: 



n 

e' 



= pj^(i-*)*a-af) i (5 . 19) 



(=1 
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e i25sziZj e i26sziZk e i2dszjZk (5 20) 

The complete adiabatic evolution is thus finally discretized in terms of the sequential action of the 
above one and two-qubit gates. 



5.2.2 Numerical results of a simulation with matrix product states 

The exact simulation of a quantum computer using matrix product states is then completely defined. 
As we said before, we have chosen the local truncation scheme in order to implement our algorithm. 
It is possible to see that the total running time of the simulation algorithm scales as 0(Tnmx 3 ). This 
reasonable truncation carries, though, an inherent - but always under control - loss of norm of the 
quantum state, since the sum of the retained squared eigenvalues will not reach 1. As we shall see, 
larger ^'s allow for more faithful simulations, as expected. 

We have implemented a number of optimizations upon the above basic scheme which are worth 
mentioning. For any non-local gate there is an overhead of SWAP operations that damage the preci- 
sion of the computation. To minimize this effect, every three-qubit clause is operated as follows: we 
bring together the three qubits with SWAPs of the left and right qubits keeping the central one fixed 
and, then, we operate the two-qubit gates. Before returning the qubits to their original positions we 
check if any of them is needed in the next gate. If so, we save whatever SWAP may be compensated 
between the two gates. Ordering of gates is also used to produce a saving of ~ 2/3 of the naive 
SWAPs. Diagonalization of the relevant reduced density matrices in the allowed Hilbert space of 
minimum dimension is used as well. A further improvement is to keep a both dynamical and local x, 
so that ancillary indices at the different partitions are allowed to take independent values and grow 
up to site-dependent and time-dependent limits. This last procedure, though, has shown essentially 
no big improvement upon a naive fixed x strategy. Let us now explain in what follows the different 
results of our simulations. 



Instantaneous expected energy 

We first simulate the adiabatic algorithm with the requirement that the right solution is found for a 
typical instance of n - 30 qubits with m - 24 clauses and T = 100. Along the evolution we compute 
the expected value of the Hamiltonian of the system, which can be calculated in C?(y 3 ) time. Our 
numerical data are shown in Fis l5.7l The system remains remarkably close to the instantaneous 
ground-state all along the approximated evolution and, as we can see, the maximum absolute error 
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with respect to our best classical simulation (x = 40) comes when evolving close to the quantum 
phase transition point. We also see convergence in the error while the system approaches the critical 
point. This minimum absolute error in the ground-state energy is, when close to criticality, of the 
order of 10~ 2 - 10~ 3 , smaller than the typical value of the energy gap for 30 qubits - as hinted by 
extrapolating the data from Fig |4.7l in Chapter 4 -. A bigger^ may bring a better precision by using 
a larger, but eventually affordable, time cost in the simulation. 



0.4 




-0.15 I ' ' ' ' 1 1 1 1 1 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

s 

Figure 5.7: Computation of the absolute error, as compared to fhe;f = 40 case, of the expected value 
of the Hamiltonian (in dimensionless units) along the adiabatic evolution for a typical instance with 
30 qubits and 24 clauses for T - 100 as^ increases. Note the increasing precision with larger ;f as 
s approaches the phase transition from the left-hand-side. In the inset, the instantaneous expected 
energy is plotted (in dimensionless units). A similar behavior is also obtained for other instances, 
getting perfect solution at the end of the computation (zero energy). 

The error in the expected energy is minimized as x increases. It is noteworthy to observe how 
the error in the simulation of the adiabatic algorithm increases at the phase transition point. We have 
also numerically checked in our simulations that it is precisely at this point where each qubit makes 
a decision towards its final value in the solution. Physically, the algorithm builds entanglement up 
to the critical point where the solution is singled out and, thereon, the evolution seems to drop the 
superposition of wrong states in the register. 

Loss of norm 

We plot in Fig l5.8l the norm of the quantum state at the end of the simulation as a function of x m 
logarithmic scale, for typical instances of 14, 18, 22 and 30 qubits. The remarkable fact is that some 
observables, like the energy, appear to be very robust against this inaccuracy, while the behavior of 
this norm was already expected not to be good, since this is precisely the parameter in which we are 
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Figure 5.8: Final norm in the register as a function of x m logarithmic scale, for instances of 
14, 18,22 and 30 qubits. 

truncating with respect to the exact evolution, and furthermore its accumulation is multiplicative as 
time evolves. 



Decay of the Schmidt coefficients 

Our simulations also allow to compute the decay of the^ Schmidt coefficients , a = 1,2, ... ,x, 
at any site / and at any step of the computation. At the closest point to criticality, and for the central 
bipartition of the system, these can be approximately fitted by the law log 2 (^ ) = a + -4= + c s/a, 
with appropriate instance-dependent coefficients a, b and c. The behavior for a typical instance of 
30 qubits is shown in Fig l5.9l 

100-qubit instance 

The ultimate goal of finding the correct solution appears also to be very robust in the simulations 
we have performed. The exact probability of success can be calculated in 0(x 2 ) time as well. As 
a symbolic example, our program has solved an instance with n = 100 bits, that is, the adiabatic 
evolution algorithm has found the correct product state out of 2 100 ~ 10 30 possibilities for a hard 
instance with m = 84 clauses and T = 2000. The simulation was done with a remarkably small 
value of x — 14 <sc 2 50 = Xmax and is presented in Fig J5. 101 Notice that while the entanglement 
entropy shows fluctuations because it is directly related to the truncation parameter of the simulation, 
the probability of success follows a smooth behavior, being almost zero at the beginning of the 
evolution, and jumping directly to one precisely when close to the quantum critical point. 



n=30 + 

n=22 x 

n=18 * 

n=14 □ 
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Figure 5.9: Decay of the Schmidt coefficients for a typical instance of 30 qubits in logarithmic scale, 
with^ = 40. The behavior seems to be approximately described by a law of the kind log 2 (/li"' /2) ) = 
a + -4^ + c y/a, for appropriate coefficients a, b and c (solid line). 




Figure 5.10: Entanglement entropy of a bipartition and probability of being at the correct solution 
as a function of s for a simulation with^ = 14 of adiabatic evolution solving a hard instance of 
n = 100 bits and m = 84 clauses. 
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Figure 5.11: Mean and worst cases of the accumulated statistics up to n - 60 for T m i n (n) (in 
dimensionless units) such that an instance is solved. Averages are performed over 200 instances for 
each n, except for n = 50, 60 with respectively 199, 1 17 instances. Error bars give 95 per cent of 
confidence level in the mean. 



Time statistics 

The robustness of evolving towards the correct solution is found for any number of qubits and small 
X- To analyze further the performance of this classical simulation, we have launched a search for 
the minimal T m j n (n) that solves samples of rc-qubit hard instances in the following way: for a set of 
small values of x, we try a random instance with an initial T, for instance T - 100. If the solution 
is found, we proceed to a new instance, and if not, we restart with a slower adiabatic evolution with, 
for instance, T = 200. This slowing down of the algorithm is performed until a correct solution is 
found and the minimum successful r mi „ is stored. Our results are shown in Fig l5. 1 II The average 
over «-qubit instances of T m i„(ri) appears to grow very slowly with n, though the extreme cases 
need increasingly larger times up to n = 25. The slowing-down in the plots for a large number of 
qubits is a side-effect of the inherent difficulty to generate hard instances of Exact Cover for large n. 
We want to remind as well that finding an instance that needs a very large T m i„ is no counterproof 
for the validity of the adiabatic algorithm, as alternative interpolating paths may solve the instance 
efficiently [66]. 

Solving hard classical instances by euclidean time evolution 

Independently of the fact that our simulation describes in an approximate way the behavior of an 
adiabatic quantum algorithm, we can think of it as a plausible classical algorithm for solving hard 
instances of an NP-complete problem. In fact, if our aim is to solve instances of Exact Cover, all 
that is required is a classical algorithm to find the ground-state of the problem Hamiltonian Hp from 
Eq J5. 161 A possibility is to perform an evolution in euclidean time, that is, to simulate the evolution 
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Figure 5.12: Euclidean time evolution solving a typical instance of 26 qubits with x — 6. The 
algorithm finds the correct solution much faster than the simulations of adiabatic quantum compu- 
tation. The sudden jump in the probability of success comes again at the maximal point for the 
entanglement entropy. 



driven by the non-unitary operator 

e~ Hpt . (5.21) 

The above evolution is not physical, since it is not unitary and therefore does not correctly preserve 
the probabilities as the parameter t (the euclidean time) flows. In any case, it is easy to see that 
if we have a (possibly not normalized) quantum state such that it has a non-zero overlap with the 
ground state of Hp, the action of the operator from Eq l5.21l over the state will eventually drive 
the original state towards the only fixed point of the map at t — > oo, which is the ground state 
of Hp. In practice, the action of the above operator over an equally-weighted superposition of all 
possible computational states will drive the original state towards the ground state of Hp with very 
high probability at times bigger than the inverse of the first gap of the system. This optimization 
algorithm can be easily implemented by using the same time-evolution procedures described before 
in terms of matrix product states. Evolution in euclidean time shall not be unitary, though, but this 
particularity does not affect any of the essential features of the updating and truncation schemes 
previously explained. 

The performance of the evolution in euclidean time for solving hard instances of Exact Cover 
is remarkably good, as compared to the performance of the simulation of the adiabatic quantum al- 
gorithm. This new classical algorithm finds the correct solution to the instances much faster than our 
previous simulations of adiabatic evolution. As an example, we show in Fig l5.12l the result of a sim- 
ulation for a typical instance of 26 qubits with x = 6. The behavior of the euclidean time evolution 
algorithm resembles very much the one of the adiabatic evolution, in the sense that the probability 
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of success remains very close to zero, until some specific point in the evolution is reached, where 
it jumps to one very quickly. It is also interesting to notice that this point corresponds, once more, 
to the point of maximum entanglement in the evolution, as measured by the entanglement entropy. 
Since the ground state of Hp is non-degenerate and separable, and since we begin with an equal su- 
perposition of all the possible states of the computational basis, the entropy must begin at zero and 
eventually die in zero, so it must necessarily reach a maximum at some point along the evolution. 
Remarkably, the point of maximum entropy coincides again with the jump in the probability of suc- 
cess. Note that even though the system is not evolving close to any quantum phase transition (like 
the one of the adiabatic quantum algorithm), the behavior along the evolution is very analogous to 
the one observed in those cases (compare Fig l5.12l and Fig l5. 10b . Again, maximum entanglement 
brings the correct solution to the problem, although our algorithm is entirely classical. 

5.3 Conclusions of Chapter 5 

In this Chapter we have shown that it is possible to implement approximated classical simulations of 
quantum algorithms by the use of matrix product states with controlled accuracy. More specifically: 

• We have implemented a simulation based on matrix product states of an adiabatic quantum 
algorithm solving the NP-complete Exact Cover problem. This simulation is made precise by 
means of an optimal local truncation scheme, and provides robust results for quantities like 
the expected energy or the probability of success, with a relatively small size of the involved 
matrices. 

• We have solved a hard 100-qubit instance of Exact Cover by means of a highly-truncated 
simulation of the adiabatic evolution algorithm. This classical simulation finds the correct 
product state out of 2 ~ 10 30 possibilities by using matrices whose indices range up to 
X = 14, much smaller than the necessary 2 50 for an exact simulation. 

• We have seen that the mean time that our approximated classical simulations take to succeed 
increases slowly with the number of qubits, though not a definite scaling law can be inferred 
given the inherent difficulty to generate very hard instances of Exact Cover for a large number 
of qubits. 

• Matrix product states algorithms for dynamical evolution can also be applied for simulating 
the non-unitary evolution in euclidean time, which we have shown to be a classical optim- 
ization algorithm that solves hard instances of Exact Cover much more efficiently than the 
classical simulations of the adiabatic algorithm. 

The results presented here could be extended in several directions. For instance, it should be 
possible to study the performance of the optimal non-local truncation scheme and to compare it 
with the one we have considered here. Also, the performance of a parallelization of the numerical 
code that we have considered here could be analyzed. More generically, it should also be plausible 
to extend the rigid structure of a matrix product state to other tensor networks specifically adapted 
to the particular problem or instance in consideration, much in the same way as PEPS do in (2 + 1)- 
dimensional systems [58]. Finally, the study of the performance of all the ideas exposed here but 
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with other quantum algorithms is a direction to be considered as well. For instance, it should 
be possible to see the behavior of a classical simulation of Shor's factoring algorithm by using 
matrix product states or related techniques. As we saw in Chapter 4, Shor's algorithm is yet another 
quantum algorithm which inherently makes use of an exponentially big amount of x in the number 
of qubits. The effect of truncations in that algorithm are, though, not evident. Perhaps, a classical 
simulation of Shor's quantum algorithm using the ideas of this Chapter could be a good candidate 
for a new classical factorization algorithm. 



Chapter 6 



Majorization arrow in quantum 
algorithm design 

Finding underlying mathematical structures in efficient quantum algorithms is one of the problems 
that quantum computation deals with. The fact that there is only a short list of ideas behind quantum 
algorithm design hints how difficult it is to come up with new quantum techniques and strategies 
to efficiently solve important problems. Graver's quantum searching algorithm [9] exploits calls to 
an oracle by enhancing a particular state, actually implementing a rotation in the relevant Hilbert 
space associated to the problem. Shor's factoring quantum algorithm [8] exploits the periodicity of 
an initial quantum state using a minimum of Hadamard and controlled-phase gates at the core of 
the quantum Fourier transform. Based on more general quantum mechanical principles, the idea 
of using adiabatic evolution to carry quantum computation [16] has proven suitable for performing 
Graver's algorithm and has been numerically studied as a candidate for attacking NP-complete 
problems, as we saw in Chapters 4 and 5. Also, the so-called quantum walks in continuous time 
have proven to efficiently solve a classically hard problem [179], whereas quantum random walks 
in discrete-time have proven to bring also Graver's square-root speed-up in a problem of quantum 
search [180]. Many other quantum algorithms can be mapped to the above families, being then 
based on the same basic principles. 

Some attempts to uncover the properties of quantum algorithms have already been explored. 
One relevant instance is undoubtedly the role of entanglement [49,50, 155-159], which was already 
considered in detail in the preceding two Chapters. In fact, although entanglement is a natural 
resource to be exploited in quantum algorithm design, there are known examples of faster-than- 
classical oracle-based quantum algorithms where the quantum register remains in a product state 
between calls to the quantum oracle all along the computation, though the speed-up is only by a 
factor of two [161, 181, 182]. In this Chapter we will concentrate on quite a different proposal. The 
basic idea is that there is an underlying strong majorization behavior in some quantum algorithms 
that seems to play a role as well. 

More concretely, we study the evolution in different quantum algorithms, with respect to major- 
ization, of the probability distribution arising in the evolving quantum state from the probabilities 
of the final outcomes, as introduced in [152]. We consider several families of quantum algorithms 
based on distinct properties. As a first step, we analyze the majorization behavior of the family of 
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quantum phase-estimation algorithms, comparing their performance with respect to majorization to 
that of Graver's algorithm [152], and giving also the explicit example of a slightly different quantum 
algorithm solving a hidden affrne function problem by means of calls to an oracle [161, 181, 182]. 
We also consider here the class of adiabatic algorithms [16] by studying the behavior of the adia- 
batic algorithm implementing a quantum search [9,69,70]. Efficiency is seen to depend on the 
interpolating time path taken along the evolution [66,69,70], and we observe that optimality in 
adiabatic quantum searching appears when step-by-step majorization is present. Finally, quantum 
walks provide exponential speed-up over classical oracle-based random walks [179], and again a 
manifest strong majorization behavior is detected. Let us begin, then, by considering the way in 
which we understand majorization theory as applied to the study of quantum algorithms. 



6.1 Applying majorization theory to quantum algorithms 

The way we relate majorization theory - as defined in Appendix A - to quantum algorithms is as 
follows: let |</r (m) ) be the pure state representing the register of a quantum computer at an operating 
stage labeled by m = 1 . . . M, where M is the total number of steps in the algorithm, and let N be 
the dimension of the Hilbert space. If we denote as (lOl^j the basis in which the final measurement 
is to be performed, we can naturally associate a set of sorted probabilities pi, i = I ... N, to this 
quantum state in the following way: decompose the register state in the measurement basis such 
that 

\^) = f j af\i). (6.1) 
The probability distribution associated to this state is 

pKm) = {p {r*) } p to s | fl («) |2 = | <# (m) )|2 ? (f> 2) 

where i = 1 . . . N. This corresponds to the probabilities of all the possible outcomes if the computa- 
tion were to be stopped at stage m and a measurement were performed. A quantum algorithm will 
be said to majorize this probability distribution between steps m and m + 1 if and only if [152-154] 

p<m) < pf.m+1) (6 3) 

Similarly, a quantum algorithm will be said to reversely majorize this probability distribution between 
steps m and m + 1 if and only if 

pim+i) < ^m) _ (6 4) 

If Eq 16.31 is step-by-step verified, then there is a net flow of probability towards the value of 
highest weight, in such a way that the probability distribution will be steeper and steeper as time 
flows in the algorithm. In physical terms, this can be stated as a very particular constructive inter- 
ference behavior, namely, a constructive interference that has to step-by-step satisfy a set of N - 1 
constraints - see Appendix A - at each time step. The quantum algorithm monotonically builds up 
the solution by means of this very precise reordering of the probability distribution. 
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It is important to note that majorization is checked on a particular basis. Step-by-step majoriz- 
ation is, then, a basis-dependent concept. Nevertheless there is a preferred basis, namely, the basis 
defined by the final measurement of the quantum register. This typically (though not necessarily 
always) corresponds to the computational basis of the quantum computer. The principle we analyze 
is rooted in the physical and practical possibility to arbitrarily stop the computation at any time 
and perform a measurement. Generically speaking, we analyze the majorization properties of the 
probability distribution of the possible outcomes of our measurement apparatus along the time-flow 
in the algorithm. 

Natural majorization 

Let us now define the concept of natural majorization for quantum algorithms. Working with the 
probability amplitudes in the basis as defined in Eq l6.ll the action of a generic unitary gate 

at step m makes the amplitudes evolve to step m + 1 in the following way: 

c/r +l) = f; U tJ af , (6.5) 

7=1 

where U{j are the matrix elements in the chosen basis of the unitary evolution operator. By inverting 
this evolution, we can write 

af = £c«af +1) , (6.6) 

7=1 

where Cjj are the matrix elements of the inverse unitary evolution, which is of course unitary as 
well. Taking the square-modulus we find 

N 

| a (m)|2 _ ^ \Cij\ 2 \af +l) \ 2 + interference terms . (6.7) 

7=1 

Should the interference terms disappear, majorization would be verified in a "natural" way between 
steps m and m+ 1 since the initial probability distribution could be obtained from the final one just by 
the action of a doubly stochastic matrix with entries \Cu\ . We shall refer to this property as "natural 
majorization": majorization which naturally emerges from the unitary evolution due to the lack of 
interference terms when making the square-modulus of the probability amplitudes. Similarly, we 
can define the concept of "natural reverse majorization", which follows in a straightforward way: 
there will be "natural reverse majorization" between steps m and m + 1 if and only if there is "natural 
majorization" between steps m + 1 and m. As we shall see, this very specific kind of majorization 
shall appear in some of our forthcoming calculations. 

6.2 Majorization in quantum phase-estimation algorithms 

Quantum phase-estimation algorithms [2, 8, 161, 181-183] are a good example of a wide class of 
quantum algorithms to begin our study. Their key ingredients are the use of the quantum Fourier 
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transform operator and the promise of a specific structure of the initial state. In [152], it has been 
numerically checked that the canonical form of the quantum Fourier transform majorizes step-by- 
step the probability distribution attached to the computational basis. Here we analytically address 
this problem and provide a proof of how the notion of majorization formulated in [152] explicitly 
operates in the special case of quantum phase-estimation algorithms. To be more specific, our 
purpose now is to present a detailed proof of the following proposition: majorization works step-by- 
step in the quantum Fourier transform of quantum phase-estimation algorithms. The whole property 
is based on the idea that Hadamard operators act by majorizing the probability distribution given 
the symmetry of the quantum state, and such a symmetry is partially preserved under the action of 
both Hadamard and controlled-phase gates [152]. 

6.2.1 The quantum phase-estimation algorithm 

Quantum phase-estimation algorithms were originally introduced by Kitaev in [183], and the basic 
problem that they aim to solve can be stated as follows. Given a unitary operator U and one of its 
eigenvectors \<p), estimate the phase of the corresponding eigenvalue U\<p) = e~ 2m ^\^>), <p e [0, 1) 
up to n bits of accuracy. An efficient solution was found in [161] and can be summarized in the 
following steps, as represented by the quantum circuit of Fis l6.ll 

(i) Prepare the pure state \\/r") = |0, 0, . . . , 0)\(p), where |0, 0, . . . , 0) is called the source register 
state of n qubits and \<p) is the target state where we have stored the given eigenvector of the unitary 
operator U. 



(ii) Apply Hadamard operators 



(6.8) 



over all the qubits i in the source state, i— 1, 2, . . . , n. 

(iii) Apply bit-wise controlled U J gates over the target state as shown in the Fig l6.ll where each 
U-i gate corresponds to the application of j times the proposed Legate with j — 0, 1 . . .n — 1. 
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Figure 6. 1 : Quantum circuit for the quantum phase-estimation algorithm. 
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Figure 6.2: Canonical decomposition of the quantum Fourier transform operator. By Uj we denote 
the unitary gate |0)(0| + e 2m/2J \l)(l\, to be controlled j - I qubits below. 



(iv) Apply the quantum Fourier transform operator 

2"-l 

2" 



QFT\q) = ^ 2 Y j e 2niqq ' ir W) 



(6-9) 



q'=0 



over the source register state. 



(v) Make a measurement of the source state of the system. This provides with high probability 
the corresponding eigenvalue of U with the required precision. 



6.2.2 Analytical results 

Let us now go through the steps of the algorithm focusing on how the majorization of the considered 
set of probabilities of the computational states evolve. The application of the Hadamard gates in step 

(ii) to the initial state produces a lowest element of majorization by means of step-by-step reverse 
majorization, 

2"-l 

1^")) = 2 -»/2 £ \x)\<p) , (6.10) 

x=Q 

yielding the probability distribution = 2~ n Vx. The outcome of the controlled £/ ; gates in step 

(iii) is the product state 

\i)/ m ) - 2~ n/2 (|0> + e" 2?n ' 2 "^|l>) • • • (|0> + e" 2,r! ' 2 >|l>) \<f>) 

= 2~ nl2 ZTJo 1 e- 2ni **\x)\<p) . (6. 1 1) 

Since the action of these gates adds only local phases in the computational basis, the uniform distri- 
bution for the probabilities is maintained (p^ - 2~ n Vx). 
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Verifying majorization for the global action of the quantum Fourier transform is simple. After 
step (iv) the quantum state becomes 

2"-l 

|^»')> = T n ^ e- 2nix ^- yl2n) \y)\(f>) . (6.12) 

x,y=0 



We then have the probability distribution 



(iv) 
Py = 



2"-l 2 



2-Z 



e 

x=Q 



-2xix(<p-y/2 n ) 



Vy . (6.13) 



Global majorization between steps (ii) and (iv) holds [152]. The remaining step (v) corresponds to 
a measurement whose output is controlled with the probability distribution p y lv \ 

While global majorization of the probability distribution is somehow straightforward to see, 
step-by-step majorization is less obvious. To this aim, the mathematical result that we shall prove 
reads as follows: the quantum Fourier transform majorizes step-by-step the probability distribution 
calculated in the computational basis as used in the quantum phase-estimation algorithm. This 
fact is seen to emerge from two important properties. It is, first, essential that the initial state 
entering the quantum Fourier transform has a certain symmetry to be discussed. Second, the order 
of the action of Hadamard and controlled-phase gates maintains as much of this symmetry as to be 
used by the rest of the algorithm. To be precise, Hadamard gates take the role of majorizing the 
probability distribution as long as some relative phases are properly protected. Controlled-phase 
transformations do preserve such a symmetry, as we shall see. 

The above property arises in three steps: the first one consists on a majorization lemma, the 
second one is a lemma concerning the preservation of phases, and finally the third one is the analysis 
of the controlled-phase operators in the quantum Fourier transform. As hinted above, we shall 
observe that the only relevant operators for the majorization procedure are the Hadamard gates 
acting over the different qubits, while controlled-phase operators, though providing entanglement, 
turn out to be immaterial for majorization purposes. 



A majorization lemma 

Let us first introduce the concept of "H(j')-pair", central to this discussion. Consider a Hadamard 
gate Ii)!} acting on qubit j of the quantum register. In general, the quantum register would cor- 
respond to a superposition of states. This superposition can be organized in pairs, each pair being 
characterized by the fact that the Hadamard operation on qubit j will mix the two states in the pair. 
Let us illustrate this concept with the example of a general quantum state of two qubits: 

|(A> =a|00>+/3|01> + y|10> + £|ll> 

= (a|OO> + 7|lO» + G0|Ol> + <5|ll» 
v v ' v v ' 

H(0)-pair H(0)-pair 

= (or|00> +/3|01» + (y|10> + 6\U)) . (6.14) 

H(l)-pair H(l)-pair 
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The second line corresponds to organizing the state as //(O)-pairs, because each pair differs only on 
the Oth qubit value. The third line, instead, organizes the state on //(l)-pairs, since each pair differs 
only on the first qubit value. We now formulate the following lemma: 

Lemma 6.1: Let denote a pure quantum state ofn qubits, with the property that the probab- 
ility amplitudes of the computational H(j)-pairs differ only by a phase for a given qubit j. Then, the 
probability distribution resulting from U^\^j) in the computational basis majorizes the one resulting 
from 

Proof: The state \if/) can always be written as: 

\iff) = ailO, 0,...,(V',...,0> + 01^10, 0,...,lA...,0> 
+ • • • + a 2 „-i|l, 1, ... ,0>, .... 1> + a 2n -ie i5 2^ |1, 1, . . . , V, ..A). 



(6.15) 



The above expression makes it explicit that the amplitudes for every pair of states that can be mixed 
by a Hadamard transformation on the qubit j only differ by a phase. The Hadamard gate JJ^ will 
mix all these pairs. The two states in every pair are equal in all their qubits except for the jth one. 
After the application of the we have 



C/^'V) = 2- 1/2 («i (l + e' 6 ') |0, 0, ... , j , . . . , 0) + a x (l - ) |0, V, . . . , 0> 

+ . . . + (l + e i6 ^) |1, 1, ... ,0>, .... 1> + a 2 „-, (l - e iS ^) |1, 1, .... 1', ... , 1>) 

We have to find a set of probabilities pk and permutation matrices P\ such that 

|ai| 2 (l+cos(5i)) 
|ai| 2 (l-cos(5i)) 

|a2"-i| 2 (l +cos(52"-i)) 
|a 2 «-il 2 (l - cos (6 2 n-i))) 

and the unique solution to this probabilistic mixture is 



(6.16) 
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(6.17) 



Pi = 



Pi = P2 = \ 

^ (0 1 

1 

Pi = 



U 



1 

1 oj 



(6.18) 



The permutation matrix P\ is nothing but the identity matrix and P 2 is a permutation of the probab- 
ilities of each pair which has undergone Hadamard mixing. This completes the proof of the lemma. 
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The lemma we have just presented states that Hadamard transformations do order the probability 
distribution when the input state has a special structure, namely, those amplitudes to be mixed only 
differ by phases. This is the key element pervading in the quantum phase-estimation algorithm: 
Hadamard transformations and controlled-phase transformations carefully preserve such a structure 
when needed, as we shall now see. 

A phase-preservation lemma 

Let us now prove the following lemma: 

Lemma 6.2: Consider the Hadamard gate acting on qubit j, and the quantum state \^ U1) ) 
from Ea \6.11\ with the property that the probability amplitudes of the computational H(i)-pairs differ 
only by a phase which only depends on i, Vi. Then, the quantum state t/^ l^'" 1 ) is such that the 
H( i)-pairs differ only by a phase Vi ± j. 

This lemma implies that the quantum Fourier transform works in such a way that states to be 
mixed by Hadamard transformations only differ by a phase all along the computation, until the 
very moment when the Hadamard operator acts. In other words, the structure of gates respects the 
relative weights of the H(/)-pairs. 

Before proving the Lemma 6.2 let us build some intuition by considering first an example. We 
start by introducing a new notation for the phases appearing in the source quantum state of Eq l6.1 II 
to be operated by the quantum Fourier transform operator by defining f3 x = -2nx(p. Then 

2"-l 

= 2 -»/2 ^ (P*\x) . (6.19) 

x=0 

Notice that since x = YICq x & 1 , we can write 

n— 1 n— 1 

p x = ^ -Inxil^ = x m , (6.20) 

;=0 i=0 

where a,- = -2n2'(p. As an example of this notation, let us write the state \ifj (m ^) in the case of three 
qubits: 

= ^ (|000> + ^ 2 |100> + e tel |010> + e i(a2+ai) \U0)) 

+ ^(|001> + e' ff2 |101> + ^ 1 |011> + / ff2+£! ' l) |lll>)e' ff0 . (6.21) 

We have factorized the ao phase in the second line of the above equation. Alternatively, we can 
choose to factorize a\, 

\if/ m ) = ^ (|000> + ^ 2 |100> + e ia °\00l) + e i{a2+ao) \m)) 

+ ^(|010> + ^ 2 |110> + e' Q ' |011> + ^ 2+£! ' o) |lll>)e" ! ' 1 , (6.22) 
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or Q-2, 

\iff m ) = (|000> + e iai |010> + e ia °\00\) + e i{a>+ao) \OU)) 

+ ^(\\00) + e iai \llO) + e ia °\lO\) + e i{ai+ao) \n\))e ia2 . (6.23) 

On the whole, the initial state for three qubits can be factorized in these three different ways. This 
example shows that there are three different possibilities to write the quantum state by focusing on 
a particular qubit. The above property is easily extrapolated to the general case of n qubits: we can 
always write the quantum state \ifr- m >) in n different ways by factorizing a particular phase in the 
second line. 

Proof: In the general case we can factorize the aj phase so that the pure state is written as 
\iff m ) = r^j (|0, 0, . . . , V, . . . , 0> + • • • + e^i^l, 1, .... o/, ... , 1>) 

+ O ' ' • • • ' lJ > ■ ■ ■ »°> + ' ' ' + e iZk * iak \h 1, .... l y , ■ • ., 1>)^ ■ (6-24) 

Then, the action of ITfi transforms the state as follows: 

u<J M m ) = ^TTvF <> <> ; ())-••• + e^m, 1 / 1>) 

+ 2 (n + l)/2 A ' + g' £ ^ a *H, 1 1 ; 1)) ■ 

(6.25) 

The resulting state still preserves the necessary symmetry property to apply Lemma 6.2 to the rest 
of qubits i ± j. The reason is that the effect of the operator has been to split the quantum state in 
two pieces which individually retain the property that all the H(/)-pairs differ only by a phase for 
i ± j. If we now apply another Hadamard operator over a different qubit, for instance qubit j — 1 , 
each of these two quantum states splits in turn in two pieces 

(1 +e ia Q(l +e' g J- 1 ) 

2(»+2)/2 

(1 +e ia J)(l -e ia J- 1 ) 
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2(»+2)/2 

(6.26) 

where ySjf is the phase defined in Eq l6.20l for the «-bit string x = (1, 1, . . . , ;_1 , 0- 7 , . . . , 1). The re- 
gister now consists of a superposition of four quantum states, each one made of amplitudes that only 
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differ by a phase. Further application of a Hadamard gate over yet a different qubit would split each 
of the four states again in two pieces in a way that the symmetry would again be preserved within 
each piece. This splitting takes place each time a particular Hadamard acts. Thus, all Hadamard 
gates operate in turn producing majorization while not spoiling the symmetry property needed for 
the next step. This completes the proof of the phase-preserving Lemma 6.2. □ 

Analysis of the controlled-phase operators 

It is still necessary to verify that the action of controlled-phase gates does not interfere with the 
majorization action carried by the Hadamard gates. Let us concentrate on the action of U^~ l \ 
which is the first Hadamard operator applied in the canonical decomposition of the quantum Fourier 
transform. Originally we had 

W m ) = -^ji (|0, 0) + • • • + ak \0, 1, . . . , 1>) 

+ J 7 j(|l,0,...,0) + --- + e , ' 2 **- ia *|l > l,...,l>)e fa - 1 , (6.27) 

where we have taken the a n -\ phase-factor out. After the action of U^~ r) we get 

,()> + ■•• + e iZk *"- ,ak \0, 1,...,1>) 
,0) + • • • + e' z **"-> ak \l, 1, . . . , 1>) = \a) + \b) . 

(6.28) 

We repeat our previous observation that the state resulting from the action of U£~ can be con- 
sidered as the sum of two states, which we have called \a) and \b). For each of these two states the 
amplitudes of the //(z')-pairs Vi + n — 1 still differ only by a phase. 

We can now analyze the effect of the controlled-phase operators. Following the structure of 
the quantum Fourier transform operator (see Fis l6.2t we focus on what happens after applying 
a general controlled-phase operator on the (n - l)th qubit of the quantum state U£~ ^iffi®) (the 
following procedure is easily extrapolated to the controlled-phase operators acting over the rest of 
the qubits). If the control qubit is the /th one, / + n - 1, then the operator will only add phases over 
those computational states from Eq l6.28l such that both the (n - l)th and the /th qubits are equal to 1, 
so we see that it will only act on the \b) state. Let us write \b) by factorizing the /th phase as follows: 

, l , . . . , 0> + • • • + i ak 1 1 , 1, . . . , Of, . . . , 1 >) 
, l l , . . . , 0) + • • • + e ! '2w,*-i ff *|l, 1 I 1 , . . . , 1)) e ia ' . 

(6.29) 

It is now clear that the action of the controlled-phase gate only adds a global phase in the second 
piece of \b), which can always be absorbed by means of a convenient redefinition of the phase 
a/. Hence we see that no relevant change is made in the quantum state concerning majorization, 



2(«+D/2 V u ' u '' 
2(»+i)/2 V 
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(l-e""- 1 ) 



+ 
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because the amplitudes of the computational //(/)-pairs V/ + n - 1 still differ only by a single phase 
which only depends on i. The action of controlled-phase operators only amounts to a redefinition 
of phases, which does not affect the necessary property for the Lemma 6.1 to hold. We see that the 
needed phase redefinition can be easily made each time one of these operators acts over a particular 
qubit. 

From all the above considerations and lemmas, it immediately follows that the quantum Four- 
ier transform operator majorizes step-by-step the probability distribution in phase-estimation al- 
gorithms, as we wished to show. We wish to emphasize the fact that controlled-phase operators play 
no role on majorization, though they provide entanglement. On the contrary, local Hadamard operat- 
ors act exactly in the complementary way, providing majorization without providing entanglement. 
We also note that the majorization arrow in the quantum algorithm is based on two ingredients. On 
the one hand we have the special properties of the quantum state, and on the other hand we have 
the structure of the quantum Fourier transform. A quantum Fourier transform acting on an arbitrary 
state would fail to obey majorization. 

One may be tempted to say at this point that Shor's quantum factoring algorithm [8] obeys 
a majorization arrow, since it can be completely understood in terms of a certain quantum-phase 
estimation algorithm, as we already saw in Chapter 4 (see Fig |4.2t . Notice, though, that there is 
a subtle but key difference between the quantum phase-estimation procedure explained here and 
the one being used in Shor's algorithm, namely, the target register in Shor's algorithm is not in a 
particular eigenstate of the unitary operator of Eq l4.4l but in a given superposition of all of them. 
This difference makes step-by-step majorization in Shor's quantum factoring algorithm fail. To see 
how this actually happens, let us remind that in Shor's quantum factoring algorithm the source state 
to be processed by the QFT operator is not the one from Eq l6.HI but the state 



for a particular / = 0, 1, . . . , r — 1 (or a superposition of all of them according to Eq l4.8b . where 
r is the period of the modular exponentiation function f(x) = a* mod N, with a randomly chosen 
a e [1,./V], N being the number to be factorized. The number of qubits n of the source register 
is chosen such that 2" € [N ,2N ]. The non-trivial instances of Shor's algorithm come whenever 
r is both even and O(N), as we saw in Chapter 4. We notice that whenever r is even, then ir + I 
is either even if I is even, or odd if / is odd, Vi. Therefore, the single bit that determines the 
parity of ir + I will always be either or 1, which implies that the corresponding qubit will always 
be either |0> or |1) in all the states of the superposition from Eq l6.30l It is clear, then, that the 
action of a Hadamard gate on that specific qubit does not majorize the probability distribution of 
the final outcomes. Even in the case of removing that qubit from the register, there typically are 
other qubits in the quantum state from Eq l6.30l that have the same value in all the states of the 
superposition, as happens already in the simple case r = 4, and about which we can not have any a 
priori information. The whole computation must be then carried without the possibility of removing 
these qubits, whose evolution breaks step-by-step majorization. Nevertheless, majorization seems 
to be working locally in the neighborhood of the final peaks of the distribution rather than globally 
on the whole set of probabilities. As a matter of fact, it is also true that our derivations rely very 
much on the specific decomposition of the quantum Fourier transform in terms of individual gates. 




i=0 



(6.30) 
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The underlying quantum circuit is not unique and majorization may not be present if alternative 
decompositions are considered. 

6.2.3 Natural majorization and comparison with quantum searching 

We now turn to investigate further the way majorization has emerged in the quantum phase-estimation 
algorithm as compared to majorization in other quantum algorithms, such as Graver's searching al- 
gorithm [9,152]. 

For a search in an unstructured database of a particular item, the best known classical algorithm 
takes asymptotically 0(2") steps in succeeding (where 2" = N is the number of entries). How- 
ever, and as we already said in Chapter 4, Graver was able to discover a quantum mechanical al- 
gorithm that implements a quadratic speed-up as compared to the best classical one, that is, Graver's 
quantum algorithm makes use of 0(2'^ 2 ) steps. We do not enter here into precise details about the 
construction of this quantum algorithm, and will only make a few comments on the way it proceeds. 
The interested reader is addressed to [9]. 

The analysis of Graver's algorithm can be reduced to a two-dimensional Hilbert space spanned 
by the state we are searching |xo) and some orthogonal state \x^ ) [2]. The unitary evolution of the 
quantum state is given by the repeated application of a given kernel K which amounts to a rotation 



where cos (0) = 1 - 2/2". Other choices of kernels are possible but the one from the above equation 
is optimal [162, 184]. The initial state of the computation is an equal superposition of all the com- 
putational states, written as = 2~ h/2 \xq) + (1 - 2~ n ) l/2 \x^) in this two-dimensional notation. For 
a given intermediate computation step the state (a,/3) T will be transformed to (a',/3') T . If we wish 
to express the initial amplitudes in terms of the final ones, we have: 



If the interference terms were to vanish then majorization would follow in a straightforward way 
from the above relations. But it is not the case. Yet it has been proven that step-by-step majorization 
in Graver's algorithm exists [152], although the way it arises is not so directly related to the unitary 
evolution in the way suggested here. 

Let us turn back to majorization in the quantum phase-estimation algorithm and its relation to 
unitary evolution. We write a generic ra-qubit state \tf/) to be operated by a Hadamard gate acting on 
the j'th qubit as 




(6.31) 



(a\ _ ( a' cos (0) + {3' sin (0) 
[pj ~ \-a' sin (0) + /?' cos (0) 



(6.32) 



We now take the square-modulus of the amplitudes, obtaining: 



\a\ 2 = cos 2 (0) \a'\ 2 + sin 2 (0) \p'\ 2 + 2 cos (0) sin (0) <Re(a'*/?) 
\p\ 2 = sin 2 (0) \a'\ 2 + cos 2 (0) \p'\ 2 - 2 cos (0) sin (0) Re(a'*/3') . 



(6.33) 



| t A> = co|0,0,...,0^,...,0> + c y |0,0,...,l^,...,0> 
+ . . . + C2 „_ w |l, 1, ... ,0>, .... 1> + c 2 n-i\\, 1, . . . , V, . . . , 1) , (6.34) 
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where we are focusing on the coefficients of the different //(j')-pairs. Applying the Hadamard gate 
over the jth qubit we get 



£/gV> = 2~ l ' 2 (co + cj)\0, (F 0) + 2- l/2 (c - cj)\0, 0, , . . , V, . . , , 0) 

+ ■ • ■ + 2- x l 2 {c 2 ->-i-j + c 2 «_i)|l, 1, .... 0/, .... 1> 
+2- 1/2 (c 2 «-i-y - C2»_i)|l, 1, .... lA .... 1> ■ (6-35) 

For a given pair of original amplitudes c m _ 7 and c m we now find final amplitudes c' . and c' m to 
be related to the initial ones as follows: 



Cm I 2 \ C m _j ~ C n 



(6.36) 



Taking the square-modulus of the amplitudes in the above expression we have 

\Cm-j\ 2 = ^IC/ + ^J 2 + ^(C;0 



C„ 



2 m ~ J 2 

2 " ^\c' m . j \ 2 + ^\c' m \ 2 -'Re(c':.jc' m ). (6.37) 



As in the Graver's previous example, we observe that if interference terms disappeared major- 
ization would arise from this set of relations. In such a case, we would only have to choose the set 
of probabilities and permutation matrices given in Eq l6.18l to prove majorization. For those terms 
to vanish, very specific properties for the coefficients c m _; and c m must hold. It can be checked that 
the interference terms vanish if and only if 



c m = a m -je i6 "'-> , (6.38) 

where real. 

The above case is indeed the case of quantum phase-estimation algorithms. Recalling our pre- 
vious lemmas, it is possible to see that the interference terms vanish also step-by-step, and therefore 
step-by-step majorization arises as a natural consequence of the unitary evolution of the algorithm. 
Notice that the quantum state from Eq l6.11l has a very specific structure so that natural majorization 
is verified step-by-step along the evolution through the quantum Fourier transform circuit. In a way 
we can say that previous steps in the algorithm prepare the source state in this particular and unique 
form, in order to be processed by the QFT operator. 



6.2.4 The quantum hidden afflne function determination algorithm 

We now wish to see how all the above properties work in a specific example of quantum algorithm, 
namely, we study majorization in a quantum algorithm solving a particular hidden affine function 
problem [181] as a generalization of Deutsch's problem [185]. Further studies have provided a 
range of fast quantum algorithms for solving different generalizations [161, 182]. The case that we 
present here is one of the multiple variations that appear in [182], but our main results are also 
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valid for the whole set of quantum algorithms that solve similar situations. As we shall see, this 
algorithm can indeed be understood in terms of a slight variation of the general quantum phase- 
estimation algorithm previously discussed. 

Let us consider the following problem [182]: given an integer- function / : Z;v — * Zjy, 
f(x) = mx + b, where x, m, b e Zjv, find out the value of m. A classical analysis reveals that no 
information about m can be obtained with only one evaluation of the function /. Conversely, given 
the unitary operator U / acting in a reversible way such that 



U f \x)\y) = \x)\y + f(x)), 



(6.39) 



- where the sum is to be interpreted as modulus N - there is a quantum algorithm solving this 
problem with only one single query to U /. The requested quantum algorithm proceeds as follows: 
let us take N - 2 n , n being the number of qubits. Perform then the following steps: 

(i) Prepare two «-qubit registers (source and target) in the state |0, 0, . . . ,0)|i/o), where \tp\) - 
QFT \\, 1, . . . , 1), and QFT denotes the inverse quantum Fourier transform in a Hilbert space 
of dimension N. 

(ii) Apply the operator QFT over the source register. 

(iii) Apply the operator U f over the whole quantum state (source and target registers). 

(iv) Apply the operator QFT~ X over the source register. 

(v) Measure the source register and output the measured value. 

The different steps concerning this process are summarized in Fia l6.3l 



|0>®« 

l«Ai) 



(n) 



(n) 



QFT 



QFT 



-l 



Figure 6.3: Quantum circuit solving the hidden affme function problem. Both source and target 
registers are assumed to be respectively composed of n qubits. 

We now show how the proposed quantum algorithm leads to the solution of the problem. Our 
analysis raises observations concerning the way both entanglement and majorization behave along 
the evolution. 

In step (i) of the algorithm the quantum state is not entangled, since that the quantum Fourier 
transform - and its inverse - applied on a well defined state in the computational basis leads to 
a separable state (see, for example, [2]). That is, the quantum state |0, 0, . . . ,0)|tAi) is completely 
separable. In step (ii) the algorithm evolves through a quantum Fourier transform in the source 
register. This action leads to a step-by-step reverse majorization of the probability distribution of 
the possible outcomes while it does not use neither create any entanglement. Moreover, natural 
reverse majorization is at work due to the absence of interference terms. 
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Next, it is easy to verify that the quantum state 

. 2"-l 

l^i> = ^72 2>~ 2 ^ /2 "^ (640) 

is an eigenstate of the operation \y) — > \y + f(x)) with eigenvalue e 2m ^ x ^ 2 " . Thus, after the third 
step, the quantum state reads 



j 2"-l 
9«/2 



2"-l ,2 OT -fc/2" 

l<Ai>- (6.41) 



x=0 



The probability distribution of possible outcomes has not been modified, thus not affecting major- 
ization. Furthermore, the pure quantum state of the first register can be written as QFT\m) (up to 
a phase factor), so this step has not eventually created any entanglement among the qubits of the 
system right after the application of the quantum oracle. 

In step (iv) of the algorithm, the action of the operator QFT~ l over the first register leads to 
the state e 2mb ^ 2 "\m)\ifri). A subsequent measurement in the computational basis over the first re- 
gister provides the desired solution. Recalling our previous results, we see that the inverse quantum 
Fourier transform naturally majorizes step-by-step the probability distribution attached to the dif- 
ferent outputs. Notice also that the separability of the quantum state still holds step-by-step. This 
observation completes our analysis of this example. 



6.3 Majorization in adiabatic quantum searching algorithms 

Our aim now is to study the majorization behavior of quantum adiabatic algorithms, which were 
already considered in the two previous Chapters. Here, we choose to analyze a very specific in- 
stance of the quantum adiabatic algorithm, namely, we consider the quantum adiabatic algorithm 
that solves the problem of searching in an unstructured database. As we shall see, the effects of a 
change of path between the initial and the problem Hamiltonian imply also a change of behavior 
in the algorithm from the majorization's perspective. More concretely, those paths leading to op- 
timality in the quantum algorithm do lead as well to step-by-step majorization, while the converse 
is not necessarily true. We do not repeat here the details of how do adiabatic quantum algorithms 
work, since they were already explained in Chapter 4. We do, however, sketch a couple of its basic 
properties. 

The quantum adiabatic evolution method has been successfully applied to the searching problem 
[69,70, 166]. Let the initial state be - -j= ^x=i l x )» N being the number of entries of the database, 
and let the initial and problem Hamiltonian respectively be Hq - I-\if/)(if/\ and Hp = I- \xo)(xq\, \xo) 
being the marked state. The interpolating Hamiltonian H(s(t)) - (1 - s(t))Ho + s(t)Hp depends on 
a time-dependent parameter s(t) satisfying the boundary conditions s(0) = and s(T) - 1 , T being 
the computational time of the adiabatic algorithm. This scheme leads to different results depending 
on whether we apply the adiabatic condition globally (that is, in the whole time interval [0, T]) or 
locally (at each time t). In what follows, we consider these two situations without entering into 
precise details of the involved calculations. For further information, we refer the reader to [69, 70] 
and references therein. 
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6.3.1 Numerical results 

We have performed a numerical analysis of the way in which majorization appears in the quantum 
adiabatic searching algorithm. Our study can be divided into two parts, regarding whether we 
demand the adiabatic condition to be fulfilled either globally or locally along the evolution. 

Analysis of the fastest global adiabatic evolution 

Let us suppose that we demand the usual adiabatic condition given in Eq l4.15l of Chapter 4 to be 
satisfied globally in the whole interval [0, T]. This does not involve any particular restriction on the 
^-dependence of s(t), so we can choose s(t) = t/T, leading to a linear evolution of the Hamiltonian. 
Under these circumstances, it can be proven [69, 70] that the global adiabatic condition is verified 
provided that 

N 

T > - , (6.42) 
e 

€ being the probability amplitude of not being at the ground-state of Hp at time T. Hence, this 
quantum algorithm needs a computational time of 0(N) to hit the right solution with high probab- 
ility, so the global adiabatic searching does not lead to an increasing efficiency with respect to a 
classical searching. 

In what follows we call P+(f) the probability of being at the marked state at time t and similarly 
P-(t) the probability of being at one of the remaining N- 1 basis states different from the desired one 
at time t. Notice that, given the symmetry of the problem, P-(t) will exactly be the same for all those 
basis states different from the marked one all along the evolution. In order to analyze majorization, 
we recall the set of inequalities given in Eq lA.3l of Appendix A to be satisfied at each majorizing 
time step. Let us make the observation that the maximum probability at all times is indeed P+(t), 
while the other probabilities will remain smaller than this quantity all along the computation and 
equal to P-(t). It is possible to see that the whole set of N cumulants that arise from the probability 
distribution follows the same basic behavior as time flows. Because of that, we present here the 
behavior of the first two non-trivial cumulants P+(t) and P+(t) + P-{t), as the rest of them do not 
lead to different conclusions. 

We have performed exact numerical simulations of the quantum algorithm in the fastest allowed 
case saturating the bound from Eq l6.42l (T = j) and have found the time evolution for the two 
cumulants. The results for e = 0.2 and N = 32 are shown in Fis l6.4l From our numerical analysis we 
conclude that a naive adiabatic quantum searching process does not produce an optimal algorithm 
neither verifies step-by-step majorization. This property is observed as the two cumulants decrease 
in time for some time steps, since there are wiggles which indicate that the system is evolving too 
fast to remain close enough to the ground state, and thus not verifying step-by-step majorization 
along the flow in time. 

Analysis of the local adiabatic evolution 

The preceding global adiabatic method can be improved if we apply the adiabatic condition given 
in Eq l4.15l of Chapter 4 locally. That is, let us divide the interval [0, T] into many small subintervals 
and let us apply Eq l4.15l to each one of these subintervals individually. Taking the limit of the size 
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Figure 6.4: Quantum searching using global adiabatic evolution with parameters e = 0.2, N - 32 
and T = 160. 

of the subintervals going to zero, we find that the adiabatic restriction has to be fulfilled locally at 
each time t: 



<e V* , (6.43) 

where H\$ is the Hamiltonian matrix element between the ground state and the first excited state 
and g(t) is the energy gap between these two states, everything given at t. This is a less demanding 
condition than Eq J4. 151 and means that the adiabaticity condition must be satisfied at each infinites- 
imal time interval. It can be shown (see, for example, [69]) that proceeding in this way the function 
s{t) must have a precise form which is given by the relation 

t = ^ — — — (arctan( - \{2s - 1)) + arctan( - 1)) . (6.44) 

We can observe this dependence in Fis l6.51 in the case of e - 0.2 and N = 32. The local adiabatic 
process implies that the smaller the energy gap between the ground and first excited states is, the 
slower the rate at which the Hamiltonian changes. With this information it can be proven [69, 70] 
that the evolution time for the algorithm to succeed with sufficiently high probability is, in the limit 

7Y » 1, 

T = -J- VyV . (6.45) 
2e 

Hence, in the case of local adiabatic evolution the computational process takes 0( VA0 time, just 
as in Graver's quantum searching algorithm, obtaining an square-root speed-up with respect to the 
best classical searching. 
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Figure 6.5: Interpolating parameter s(t) for quantum searching using local adiabatic evolution. 



Defining P+{t) and P~{t) as before, we can again restrict ourselves to the study of the two non- 
trivial cumulants P+(t) and P+(t) + P~(t) in order to observe the evolution of majorization. We have 
numerically solved the dynamical equations for e = 0.2 and N = 32, and have found the evolution 
of the two quantities, which is given in Fig J6.6l From the numerical analysis, it follows that a local 
adiabatic searching algorithm is not only optimal in time, but also verifies step-by-step majorization. 



Analysis of slower global adiabatic evolutions 



Let us now consider global adiabatic evolutions which are not necessarily tight in time, that is, 
extremely slow time variations of the Hamiltonian, much slower than the minimum necessary for 
the adiabatic theorem to hold. In the case we are dealing with, this implies the consideration of the 
case in which T > — , that is, the adiabatic inequality from Eq l6.42l is not saturated. 

We have again performed a numerical analysis for the time evolution of the two non-trivial 
cumulants P+(t) and P+(t) + P-(f), for e = 0.2, N = 32, and T = 320 and 480 (both cases bigger 
than — = 160). The results are plotted in Fie l6.7l and Fie l6.8l From these two plots, we observe 
that a step-by-step majorization tends to appear as long as the evolution of the Hamiltonian becomes 
slower and slower. From a physical point of view, this means that the probability of "jumping" to 
the first excited state decreases as long as the evolution is performed at slower changing rates, thus 
satisfying better the assumptions of the adiabatic theorem. Step-by-step majorization may thus 
appear in global adiabatic searching processes for a slow enough evolution rate. 



6.3. Majorization in adiabatic quantum searching algorithms 
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Figure 6.6: Quantum searching using local adiabatic evolution with parameters e = 0.2, N = 32 and 
T = 44. 
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Figure 6.8: Quantum searching using global adiabatic evolution with parameters e - 0.2, N = 32, 
and T = 480. 

6.4 Majorization in a quantum walk algorithm with exponential speed- 
up 

The extension of classical random walks to the quantum world has been widely studied, yielding 
two different models of quantum random walks, namely, those which operate in discrete time by 
means of a "coin operator" [186-188] and those based on a Hamiltonian evolution in continuous 
time [179, 189, 190]. Regarding the discrete-time model of quantum random walk, two indicative 
algorithmic results have been found, namely, an exponentially fast time when crossing the hyper- 
cube with respect to the classical random walk [191] and a quantum searching algorithm achieving 
Graver's quadratic speed-up [180]. As a matter of fact, the first one of these two results does 
not provide any algorithmic speed-up, as there exists a classical algorithm that solves the hitting 
problem in the hypercube exponentially faster than the naive classical random walk, that is, in a 
time 0(poly(log 2 AO) where N is the number of nodes of the graph (see [191, 192]). Nevertheless, 
the second of these examples shows algorithmic advantage with respect to any possible classical 
strategy. The analysis of the quantum random walk searching algorithm shows that the quantum 
evolution can be understood as an (approximate) rotation of the quantum state in a two-dimensional 
Hilbert space which is exact in the limit of a very large database (see [180] for details), resembling 
the original proposal of Graver's searching algorithm which can be decomposed exactly in a two- 
dimensional Hilbert space (see Eq l6.32b . This rotational structure of the evolution implies again 
step-by-step majorization when approaching the marked state, exactly in the same way as the usual 
Graver's searching algorithm [9, 152]. 

Here we wish to restrict ourselves to the continuous-time model of quantum walk and analyze 
a proposed quantum algorithm based on a quantum walk on continuous time solving a classically 
hard problem [179]. We sketch the main ingredients of the problem setting and its efficient solution 
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Figure 6.9: A possible graph constructed from two binary trees with n = 4. 

in terms of a quantum evolution (the interested reader is addressed to [179] for specific details). For 
a more generic review on quantum walks both in discrete and continuous time, see [192]. 

6.4.1 The exponentially fast quantum walk algorithm 

The problem we wish to solve is defined by means of a graph built in the following way (see 
[179]): suppose we are given two balanced binary trees of height n with the 2" leaves of the left 
tree identified with the 2" leaves of the right tree in a simple way, as shown in Fig l6.9l A way 
of modifying such a graph is to connect the leaves by a random cycle that alternates between the 
leaves of the two trees, instead of identifying them directly. An example of such a graph is shown 
in Fig lOOl 

Suppose that the edges of such a graph are assigned a consistent coloring (that is, not two edges 
incident in the same vertex have the same color), and that the vertices are each one given a different 
name (with a 2«-bit string, so there are more possible names than the ones assigned). We now define 
a black-box that takes two inputs, a name a given as a 2n-bit string and a color c, and acts in the 
following way: if the input name a corresponds to a vertex that is incident with an edge of color c, 
then the output corresponds to the name of the vertex joined by that edge; if a is not the name of a 
vertex or a is the name of a vertex but there is no incident edge of color c, the output is the special 
2«-bit string (1,1,..., 1), which is not the name of any vertex. 

Now, the problem we wish to solve reads as follows: given a black-box for a graph such as the 
one previously described, and given the name of the IN vertex, find out the name of the OUT vertex. 

In [179] it was proven that no classical algorithm can transverse a graph such as the one in 
Fig J6. 101 in polynomial time, given such a black-box. Furthermore, an explicit construction of a 
quantum algorithm based on a continuous-time quantum walk on the graph that succeeds in finding 
the solution for this oracular problem in polynomial time was given. The quantum algorithm of 
[179] for this problem can be briefly summarized as follows: consider the (2n + 2)-dimensional 
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Figure 6.10: An alternative graph constructed from two binary trees with n = 4. Connection 
between the leaves is made through a random cycle. 



|colj> = -i= J] |a), (6.46) 

V J a e column j 



subspace spanned by the states 



where Nj = V if < j < n and Nj = 2 2n+l ~ j if n+ 1 < j < 2n+ 1. We call this subspace the "column 
subspace", and each state of the basis is an equally weighted sum of the states corresponding to the 
vertices lying on each column of the graph. We now define a Hamiltonian acting on this subspace 
by the following non-zero matrix elements: 



1 , if0</'<rc-l,rc+l<7<2n 
(col U + l)l#|col j) = (col j\H\co\ (j + 1)) = ' . J (6.47) 

1 2 1/z , if j - n . 

The action of this Hamiltonian on the graph is nothing but promoting transitions between adjoint 
vertices, so a quantum walk on the graph (on the whole Hilbert space) generated by this Hamiltonian 
is equivalent to a quantum walk on the line (on the column subspace). Because of that, from now 
on we only focus our attention on the quantum walk on the line generated by the Hamiltonian from 
Eq l6.47l Moreover, it can be proven that given the structure of the graph in the form of a black- 
box such as the one already described, our Hamiltonian can be efficiently simulated by means of a 
quantum circuit [179]. 

The quantum walk works as follows: at first the "wave packet" will be precisely localized at the 
IN vertex (the initial state will be |col 0)). Due to the unitary time evolution driven by the Hamilto- 
nian, it will initially spread out through the different vertices at the left hand side of the graph (those 
belonging to the left binary tree), but after a short time (once half the graph has been trans versed) 
it will begin to spread through the vertices on the right hand side, interfering constructively in the 
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Figure 6.11: Probability of finding the OUT node in the quantum walk algorithm, for n = 4. 

OUT vertex as the time goes on. Physically, this is nothing but a wave propagation. Should we 
wait longer, the wave packet would come back to be localized at IN vertex and the process would 
similarly be repeated again. Actually, due to the "defect" of the Hamiltonian in the central vertices, 
it can be shown that the transmission through the central columns is not perfect, but high enough 
for the OUT node to be achieved with a very high probability in small computational time. In [179] 
the authors prove that the succeeding time is polynomial in n. 



6.4.2 Numerical results 

We have numerically simulated this quantum walk for the particular case of n = 4, and have plotted 
the time evolution of the probability of success in Fig l6.11l We observe that the numerical result 
is in agreement with the prediction that the time the algorithm takes in hitting the OUT node with 
high probability seems to be, at first sight, linear with the size of the system. 

In order to analyze majorization, for the case n = 4 there are 62 cumulants that can be computed 
from a set of 10 non-trivial probabilities. This is so due to the fact that all the states of the whole 
Hilbert space belonging to the same column always share the same probability amplitude. The 
quantities to be considered are then the probabilities of being at each column state normalized by 
the number of nodes belonging to that column, that is, the probability of being in one node of each 
column. In general, there are then 2n + 2 different probabilities to be considered at each time step. 
Given only these 10 quantities, we were able to compute all of the 62 cumulants corresponding to 
all the partial sums of sorted probabilities, according to Eq lA.3l in Appendix A. In order to make the 
figures as clear as possible we have only plotted 10 of these quantities in Fig l6.12l which correspond 
to the cumulants arising from the sorted probabilities when only one node per column is considered. 
Our numerical simulations indicate that the rest of the cumulants exhibit a behavior similar to that 
of the ones appearing in Fig l6. 121 and thus bring no further insight. We have also numerically 
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Figure 6.12: Time evolution of the ten cumulants in the quantum walk algorithm when one node 
per column is considered, for n = 4. The evolution follows a majorization cycle 

simulated the algorithm in the case of a bigger graph, namely, in the case n = 10. In this case there 
are 2n + 2 = 22 different probabilities to be considered at each time step. Proceeding in the same 
way than in the case n = 4 (that is, not plotting all the cumulants, but the only the sorted sum of 
these 22 probabilities), we obtain a similar behavior as in the case for n = 4, as is shown in Fig l6.13l 

Looking at the two plots, we conclude that the continuous time quantum walk follows a step- 
by-step majorization cycle all along the computation until it reaches the OUT node. It is worth 
remarking as well that the time the algorithm spends reversely majorizing the probability distribu- 
tion is about half of the time of the whole computation. The physical reason for this behavior is 
clear, as this is the time the "wave packet" spends spreading over the binary tree on the left hand 
side, thus leading to a destructive interference part. Note that such a destructive interference indeed 
strictly follows a step-by-step reverse majorization of probabilities. Furthermore, by combining 
Fig l6.11l and Fig l6.12l we see that the raising of the probability of success is linked to a step-by- 
step majorization. Physically, this is the part in which the algorithm constructively interferes into 
the OUT node once the wave packet is approximately in the right-hand-side binary tree. We see 
that this constructive interference follows a majorization arrow. Actually, the observed majorization 
cycle is very similar to the one that we already found in the quantum phase-estimation algorithm, 
but in this case we have numerically checked that the present cycle does not seem to follow the rules 
of natural majorization. Complementarily, we have also observed that the probability amplitudes 
follow the rule that those belonging to even columns are real, while those belonging to odd columns 
are imaginary. 

The quantum random walk heavily exploits the column structure of the problem. The register 
works on a superposition of columns, that is of states belonging to the same column with equal 
weight. It is then natural to ask whether a step-by-step majorization cycle operates also at the level 
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Figure 6. 13: Time evolution of the 22 cumulants in the quantum walk algorithm when one node per 
column is considered, for n = 10. The cumulants tend to collapse in the plot given the size of the 
graph. The evolution follows a majorization cycle. 

of columns. The idea behind this analysis corresponds to accept that the final measurement will filter 
each one of the columns as a whole. The result of the measurement would correspond to determining 
a particular column. The point here is to find to what extent the success of finding the OUT state 
is related to the column structure of the algorithm. We have numerically considered the column 
amplitudes for n = 4 and n = 10 with a total of 9 and 21 cumulants to be calculated respectively 
from the sorted probabilities at each time step of being at each column of the graph. In Fig J6. 141 
and Fig l6. 151 we plot our results, which show that there does not exist a majorization cycle when 
the final measurement is carried on columns. The conclusion is that deterministic quantum walks 
cleverly exploit the column subspace structure of the problem to achieve step-by-step majorization 
on the individual states, but not on the individual columns. 

6.5 Conclusions of Chapter 6 

We have seen in this Chapter that majorization seems to appear in the fauna of quantum algorithms 
in a very specific way, namely, in such a way that some instances of efficient quantum algorithms 
seem to step-by-step majorize the probability distribution of the final outcomes all along the flow in 
time. In order to be precise: 

• We have proven that the quantum Fourier transform in quantum phase-estimation algorithms 
majorizes step-by-step the probability distribution of the final outcomes. This step-by-step 
majorization is seen to appear in a natural way from the absence of some interference terms 
in the unitary evolution, in contrast with what is found for Graver's quantum searching al- 
gorithm. The example of a quantum algorithm solving a hidden affine function problem also 
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Figure 6.15: Time evolution of the 21 cumulants in the quantum walk algorithm when the column 
measurement is considered, for n - 10. No majorization cycle is present. 
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shows the same basic features than the quantum-phase estimation algorithm, which can be 
understood in terms of a majorization cycle along the complete time-evolution. However, 
Shor's quantum factoring algorithm, though being based on a variant of the quantum phase- 
estimation algorithm, does not globally obey step-by-step majorization on the whole set of 
relevant probabilities. 

• We have seen that step-by-step majorization in adiabatic quantum searching algorithms is 
heavily attached to the optimality of the interpolating path. Those paths which do not pro- 
duce an optimal quantum search are seen to step-by-step majorize the probability of the final 
outcomes only if the change rate of the Hamiltonian is extremely slow. On the contrary, the 
optimal path producing a square-root speed-up directly obeys step-by-step majorization. 

• We have observed that there is a majorization cycle of the probabilities of the final outcomes 
in an exponentially fast quantum walk algorithm solving a classically hard problem defined in 
terms of a non-trivial graph. This majorization cycle does not appear if alternative collective 
measurements are considered. 

Our conclusion is that some broad families of quantum algorithms seem to have an underlying 
majorization structure in the way they proceed in order to get the desired solution to the problem 
that they deal with. This behavior is somehow similar to the one of greedy algorithms in classical 
computation, which always evolve such that the probability of the "winner" increases at each time 
step. Majorization is, though, a far more severe condition, since it not only involves constraints 
on one single and specific probability, but on the complete probability distribution. In some sense, 
majorization seems to be a plausible candidate to look at in order to have a good understanding of 
the performance of a quantum algorithm, together with entanglement. How these two quantities 
behave along the computational evolution of a given quantum algorithm may already provide a lot 
of information about its performance. 



Chapter 7 

General conclusions and outlook 



The work presented in this thesis tries to bring together different fields of physics. We used tools 
from quantum information science to analyze problems in quantum field theory and condensed- 
matter physics in Chapters 1 and 3. Conformal field theory can in turn be useful to analyze prob- 
lems in quantum information science, as we saw in Chapters 1 and 2. Moreover, quantum phase 
transitions and quantum algorithms are seen to be very much related, as we have seen in Chapters 
4 and 5. Furthermore, ideas related to the performance of some quantum algorithms were shown in 
Chapter 6 by using majorization theory. All in all, we have seen that the fields of quantum informa- 
tion science, condensed-matter physics, and quantum field theory have very much in common, and 
that their multidisciplinary intersection is useful. 

Let us consider several future directions. First, the use of majorization theory and conformal 
field theory together with related techniques applied to a comprehension of both the irreversibility 
of renormalization group flows and the behavior of the single-copy entanglement in more than 
(1 + 1) dimensions is something that remains to be done. Also, it is still a theoretical challenge to 
know whether adiabatic quantum algorithms can solve NP-complete problems in polynomial time 
or not, which in the end amounts to ask about the possibilities of quantum computation to solve 
the celebrated P^NP conjecture. Further analysis of adiabatic quantum algorithms could be done, 
for instance, by means of a parallelization of the local truncation scheme that we used in Chapter 
5, or by means of non-local truncation schemes, adapted valence-bond ansatzs for the ground state 
wavefunction, or other related techniques. Indeed, classical numerical simulations using the ideas 
from Chapter 5 of some other quantum algorithms, like Shor's factoring quantum algorithm, could 
bring further insight both for the quantum algorithm and for the classical simulation technique itself. 
The big problem in quantum computation remains to be, yet, the design of new, useful and efficient 
quantum algorithms. Furthermore, from the many-body physics point of view, the challenge now 
is to perform reliable and accurate classical simulations of the properties of (2 + l)-dimensional 
quantum many-body systems, for which new numerical techniques are beginning to be discovered. 
However, the tools developed so far do not apply to the study of critical fermionic systems in more 
than one spatial dimension, since some of these systems break the entropic area-law scaling [109- 
111]. A better understanding of these models, both from a theoretical and numerical point of view, 
together with a plausible numerical ansatz for their ground state wave function, remains as an open 
problem. 
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Appendix A 

Majorization 



Majorization theory deals with the notion of relative order of probability distributions. It was origin- 
ally introduced within the fields of mathematical statistics and economics [1 1-14], and its basic idea 
relies on the comparison of two given probability distributions by means of a set of order relations 
to be satisfied by their components. 

We now precisely define the notion of majorization [14]. Let x, y e R +N be two normalized 
probability vectors, Y^=\ x i = Z^iJi = 1- We say that distribution y majorizes distribution x, 
written x < y, if and only if there exist a set of permutation matrices and probabilities > 0, 
YikPk - 1» sucn that 



Since, from the previous definition, * can be obtained by means of a probabilistic combination of 
permutations of y, we get the intuitive notion that probability distribution * is more disordered than 
probability distribution y. This defines a partial order in the space of probability distributions. 

There are two alternative equivalent definitions of majorization which turn out to be useful. The 
first one reads as follows. We say that a given N x N matrix D is doubly stochastic if it has non- 
negative entries and each row and column adds up to 1. Then, y majorizes x if and only if there is a 
doubly stochastic matrix D such that 



Notice that in Eq lA.ll Y>k PkPk = D defines a doubly stochastic matrix, that is, D has nonnegat- 
ive entries and each row and column adds up to unity, thus satisfying Eq lA.2l 

The third equivalent definition of majorization can be stated in terms of a set of inequalities 
between partial sums of the two distributions. Consider the components of the two probability 
vectors sorted in decreasing order. Then, x < yif and only if 




(A.l) 



k 



(A.l) 



k 



k 





(A3) 



7=1 



7=1 



All along this thesis, we refer to these partial sums of sorted probabilities as cumulants. 
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A powerful relation between majorization and any convex function / over the set of probability 
vectors states that 

2<?=>f(2)<ftf). (A.4) 

From this relation it follows that the Shannon entropy H(z) = - Tif=i Zi 1°§2 z i °f a probability 
distribution zeR* satisfies H(x) > H (y) whenever x < y. Majorization is, therefore, a stronger 
notion of order for probability distributions that the one imposed by the entropy H(z). 

The connection between majorization and quantum mechanics can be established whenever a 
probability distribution appears. For instance, one could be interested in the majorization properties 
of the probability distribution arising from the spectrum of some given reduced density matrix, as 
happens often in the field of quantum information science. For two reduced density operators p and 
cr with spectrums p and cr, we say that p < cr if and only if p < cr. This extends the notion of 
majorization to positive semi-definite operators by considering their normalized spectrum. 



Appendix B 



Some notions about conformal field 
theory 



The aim of this Appendix is to give a brief, non-technical and non-exhaustive idea about some of the 
basic concepts of conformal field theory. The interested reader is referred to the specific literature 
in the field for further details and developments (see for example [21] and references therein). 

Consider a metric g^ix) of signature (p, q) in a space of total dimension D, where x stands 
for a given point of this space in some given coordinate system. Under a change of coordinates 
x —* x', the metric transforms as g'^x') - jj^§^gap(x), where sums are to be understood on 
repeated indices from now on. The conformal group in D dimensions is, by definition, the subgroup 
of coordinate transformations that leave the metric invariant up to a local change of scale, 

8iiv(x) -> g'^yix) = Q.(x)g MV (x) , (B.l) 

where Q(x) is a local dilatation factor. It is possible to exactly characterize the form of these trans- 
formations, which are given by the Poincare group 

x — > x' = x + a 

x -> x' = Ax (AeSO(p,q)) (B.l) 

with Q(x) = 1 , the dilatations 

x -> x = Ax (B.3) 
with Q = /T 2 , and the so-called special conformal transformations 

, x + bx 2 ^ „ s 

x -» x' = - — — (B.4 

1 + 2b ■ x + b 2 x 2 

with Q(x) = (1 +2b-x+b 2 x 2 ) 2 . Conformal symmetry can then be understood as some generalization 
of scale symmetry. Those field theories defined in the continuum that are invariant under conformal 
transformations constitute the so-called conformal field theories. 

Conformal symmetry is especially powerful in the case of 2 dimensions, typically denoted as 
(1 + 1), in the case of having one temporal and one spatial dimension. Given the coordinates of the 
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plane x 1 and x 2 , and defining new complex coordinates z, = x 1 + ix 2 and I = x 1 - ix 2 (respectively 
called holomorphic and antiholomorphic coordinates), conformal transformations in 2 dimensions 
coincide with the set of analytic coordinate transformations in the plane 

z -» m 

z -» fiz), (B.5) 

/ and / being analytic complex functions. Typically, it is useful to work with z and I treated as 
independent variables, so that the physical condition z, = z* is left to be imposed at our convenience. 
The fact that conformal transformations in the plane precisely coincide with the group of analytic 
coordinate transformations is very notorious, since the number of generators of the conformal group 
in 2 dimensions is then infinite, which only happens for this number of dimensions. The behavior of 
conformally-invariant field theories in 2 dimensions is, then, heavily constrained by the symmetry. 

In order to be more specific, assume that we are given a conformally-invariant quantum field 
theory in D = 2. Those operator fields <I>(z, z) that transform under conformal transformations like 

df\' l (df) h 



° (Z ' !H i) 1«) $(/(z),/ ® )i (R6) 

with positive real h and h, are called primary fields of conformal weight (h, h). Conformal symmetry 
imposes that the two-point correlation function of two primary fields (®i(zi,Zi)Q>2(z2,Z2)) must De 

<Ol(Zl,Zl)«2(Z2.Z2)> = -z^7 (B-7) 
_2ft=2ft 

if (h\, h{) = {hi,hi) and zero otherwise, where z\2 - Z\ -Zi, Z12 - Zi -Zi- Note that the decay of the 
correlation function in Eq lB.7l is algebraic, as is the typical situation of critical condensed-matter 
systems. This is not strange, since many critical quantum many-body systems can be understood at 
criticality as the regularization on a lattice of some given conformal field theory, as is the case, for 
example, of the critical Ising quantum spin chain [21,22]. Indeed, conformal symmetry imposes 
similar decaying laws for the two-point correlators in any number of dimensions. 

An important quantity which is to play a role is the stress-energy tensor T MV (x), which can be 
always defined for any field theory. For instance, for a free-bosonic quantum field theory defined in 
terms of a Lagrangian the stress-energy tensor reads 



0(0^0) 



where tp stands for the quantum field of the free boson. It can be seen that in two dimensions, the 
stress-energy tensor of a conformally-invariant quantum field theory has only two non-vanishing 
components, which are called T(z) and f(z). An important property of a primary field <l>(w, w) is 
that its operator product expansion with the stress-energy tensor reads 

h 1 „ 

T(z)®(w, w) = -r<I>(w, w) + -3 w <D(w, w) + ■ ■ ■ 

(z - w) z (z - w) 

h 1 

r(z)0(w, w) = -r<D(w, w) + -d^O(w, w) + ■ ■ ■ , (B.9) 

(z - w) 1 (z - w) 
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which can be understood as an alternative definition of a primary field of conformal weight (h, h). 

The stress-energy tensor is an example of a quantum field that is not primary. Computing its 
operator product expansion with itself, one gets 

T(z)T(w) = -^— + -^—T{w) + —^—d w T{w) 

(z - wp {z - w) L (z - w) 

f(z)f(w) = -^-7 + — ^—jf(w) + — !— d w T (w) , (B.10) 

(z - wp (z - w) 1 (z - w) 

which clearly differs from Eg IB .91 The above equations define the so-called holomorphic and an- 
tiholomorphic central charges c and c, which depend on the particular theory under consideration, 
much in the same way as the conformal weights (h, h) do. For example, for a free bosonic quantum 
field theory c — c — 1, whereas for a free fermionic quantum field theory c - c = 1/2. Yet, 
another property of the stress-energy tensor for conformally-invariant quantum field theories in 2 
dimensions is that it is possible to expand it in terms of modes as follows: 

T{z) = J]z~ n - 2 L n 

m = j]r i - 2 L„, (B.n) 

neZ 

where the operators L n and L n satisfy the commutation relations 

c 3 

[L„, L m ] = in - m)L n+m + —{n - n)6 n+m fl 

- - ~ ^3 

[L n , L m ] =(n- m)L n+m + —{n - n)5 n+m fi 

[L n ,L m ] = 0. (B.12) 

The above equations define two copies of an algebra which is called the Virasoro algebra. Every 
conformally-invariant quantum field theory determines a representation of this algebra, with some 
c and c. 

The construction of the Hilbert space for a conformal field theory in 2 dimensions is very much 
related to the above operator algebra. Given a vacuum |H) which is assumed to exist by hypothesis, 
the state 

\h,h) = <D(0,0)|Q> (B.13) 
created by a primary field 3>(z, z) of conformal weight (h, h) satisfies 

L \h,h) = h\h,h) 
L Q \h,h) = h\h,h) 

L n \h,h) = L m \h,h) = Vrc,m>0. (B.14) 
Any state satisfying the above relations is called a highest-weight state. States of the form 

L—niL—rf 2 • ■ • L- n jL- mi L- m2 ' ■ ' L- mk \h, fi) (B.15) 
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are called descendant states, and are also eigenstates of Lq and Lq with eigenvalues h+n\ +ri2+- ■ -+tij 
and h + m\ + ni2 + • • • + respectively. The full tower of eigenstates of Lq and Lq constructed in 
this way is known as the Verma module. Therefore, the Hilbert space of a conformally-invariant 
quantum field theory in 2 dimensions decomposes as the direct sum of Verma modules, the number 
of which depends only on the number of primary fields appearing in the theory. 



Appendix C 



Some notions about classical complexity 
theory 

In this Appendix our aim is to give some very basic notions and non-technical background on 
classical complexity theory. Excellent textbooks on this topic are those of Garey and Johnson [193] 
and Papadimitriou [194]. A review on complexity theory, with extensions to quantum complexity 
theory, is given by Aharonov and Naveh in [151]. 
Let us begin with the following definition: 

Definition C.l: An alphabet £ is a set of symbols. 

We did not define the concept of symbol since we believe its meaning to be clear from the 
context. Examples of alphabets are Ei = {a, b, . . . ,z}, E2 = {a, ft, ■ ■ ■ , a)}, and £3 = {0,1}. The 
alphabet S3, with only two symbols, is usually referred to as the binary alphabet. 

Definition C.2: A language L over an alphabet £ is a set of strings of symbols from £. 

For instance, L\ = {jack, sam,daniel, teak} is a language over the alphabet Ei, and L2 = 
{010, 00010, 1001} is a language over the binary alphabet E3. 

Definition C.3: A decision problem is a problem for which the answer belongs to a binary 
alphabet. 

This is the kind of "yes" or "no" problems. That is, questions of the type "will the universe 
expand forever?", or "do I prefer chocolate or lemon ice-creams?", but also questions like "is the 
number 1761935875391 the product of two or more primes?". An important part of the theory of 
computational complexity is built in terms of decision problems. More concretely, one has to decide 
whether a given string of symbols from an alphabet, called instance, belongs to a certain language 
or not. From now on we shall always restrict ourselves to the binary alphabet, whose symbols are 
called bits. 

Languages are classified in terms of complexity classes, according to different criteria. We now 
define a complexity class that plays a major role in complexity theory: 
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Definition C.4: P is the class of languages L for which a deterministic Turing machine can 
decide in a time 0(poly(|x|) if an instance x belongs to L or not, \x\ being the number of bits of x. 

In the above definition, we understand that a deterministic Turing machine is our classical model 
of computation. Usually, it is said that languages LeP can be decided in polynomial time by a de- 
terministic Turing machine. Intuitively, we understand that a language L belongs to the complexity 
class P if there is an efficient classical algorithm that allows to deterministically decide whether a 
given instance x belongs to L or not, where by the term "efficient" we mean "polynomial in the size 
of the instance". Let us now define another important complexity class: 

Definition C.5: NP is the class of languages Lfor which there exists a deterministic polynomial- 
time verifier V such that 

• VxeL, there is a y such that \y\ = poly(|x|) and V(x,y) = 1, and 

• Vx g L and Vy such that \y\ - poly(|x|), V(x,y) = 0. 

Usually y is referred to as the witness or certificate. Both the witness y and the verifier V help 
in deciding whether the instance x belongs to L or not. Let us clarify Definition C.5 by means of an 
example: let L - COMPOSITE be the language of numbers that can be decomposed as a product 
of two or more primes. Let x - 161 be an instance of the decision problem "does x belong to 
COMPOSITE?". A possible witness y can be given by the two prime numbers 7 and 23, and the 
verifier V can be a classical deterministic algorithm that performs the following check: 7x23 = 161. 
Notice then that if the instance 161 belongs to COMPOSITE there is a witness 7, 23 such that the 
verifier can check that the instance belongs to the language. On the contrary, if we are given an 
instance that does not belong to COMPOSITE (for instance, x - 17), then there is no witness y such 
that our verifier can check that 17 is a product of two or more primes. In a way, the witness has 
to be thought of as the "proposal of solution", and the verifier has to be considered as a classical 
algorithm that allows to deterministically and efficiently check whether the proposed solution to 
the specific instance is correct or not. This example shows that COMPOSITE e NP, which in less 
mathematical words is commonly referred to as "the problem of deciding whether a given number 
is the product of two or more primes is NP". 

Given the Definition C.5 of the NP complexity class, we can now define the following: 

Definition C.6: NP-hard is the class of languages L such that the problem of deciding whether 
an instance x' belongs or not to a language L' e NP can be efficiently reduced to the problem of 
deciding whether an instance x belongs or not to L, Vx' and V e NP. 

In plain words, a problem is said to be NP-hard if all the instances of all the NP problems can be 
efficiently mapped to specific instances of the NP-hard problem. Therefore, if a language L e NP- 
hard can be decided by some deterministic classical algorithm, the same procedure can essentially 
be applied to decide all the languages in the complexity class NP, and "solve all the NP problems". 

Let us now define the important concept of NP-complete: 

Definition C.7: NP-complete is the class of languages L such that L e NP-hard and L e NP. 
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According to Definition C.7, NP-complete languages are those languages in NP such that being 
able to decide about one of them implies being able to decide about the whole complexity class 
NP. An important example of an NP-complete language is 3-SAT. A possible instance of the 3-SAT 
decision problem is a boolean formula in conjunctive normal form over n bits <p(x\, x%, . . . , x n ) = 
C\ A C2 A • • • A C m , where X{, i - 1, 2, . . . , n, denotes the value of the bits, and Cj, j - 1,2,..., m, 
are the so-called clauses. Each clause Cj is built in the way Cj = (xjj V Xj^. V x/,3), where xy ;Q . is 
a literal for bit a of clause j, which can be any of the n bit variables or its negation. The decision 
problem is properly defined by the following question: "given an instance <f> is there a string of 
n bits (yi,y>2, ■ ■ ■ ,y n ) such that (p(y\,y2, ■ ■ ■ ,y n ) - 1?", or equivalently, "is there a string of n bits 
(y\,y2, ■ ■ ■ ,y n ) such that all the m clauses are satisfied?". 

The proof of the NP-completeness of 3-SAT is one of the most relevant results in the field of 
complexity theory, and is due to the original work of Cook [72]. That proof opened the door to the 
discovery of many other NP-complete languages and, today, NP-complete languages (or problems) 
appear in many different fields of mathematics, physics, and science in general. Their relevance 
comes in part from the fact that they are at the heart of one of the most celebrated open questions in 
mathematics, which reads as follows: 

Problem C.l: /sP# NP ? 

To determine the answer to the above question, it would be sufficient to prove that it is possible 
to deterministically decide some NP-complete language efficiently, and then P = NP, or on the 
contrary to prove that it is impossible to efficiently and deterministically decide an NP-complete 
language, and therefore P + NP. While the most accepted opinion is that P + NP, it has been so far 
impossible to produce a precise and mathematical proof of this, neither of the opposite statement 
P = NP. Indeed, Problem C.l remains today as probably the most challenging open problem in 
computer science [194]. 

Let us mention as well that the deterministic complexity classes P, NP, NP-hard and NP-complete 
can be further generalized if we consider classical probabilistic models of computation, the equival- 
ent probabilistic complexity classes being called BPP, MA, MA-hard and MA-complete. Further- 
more, if the underlying computational model is a quantum computer, the corresponding generalized 
quantum complexity classes are called BQP, QMA, QMA-hard and QMA-complete. The study 
of these classes is beyond the scope of this Appendix, and we refer the reader to [151] and refer- 
ences therein for further details on quantum complexity theory and its consequences for quantum 
computation. 
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